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1.0 


BACKGROUND 


As detailed in the original Statement of Work, the objective 
of Phase II of this research effort was to develop a general 
framework for rocket engine performance prediction that 
integrates physical principles, a rigorous mathematical 
formalism, component level test data, system level test data, and 
theory-observation reconciliation. Specific Phase II development 
tasks are defined as follows: 

Task 1. Identify device modules required for rocket engine 

analysis. Identify device specific forms of fundamental 
physical principles. Identify device specific forms of 
generally accepted engineering approximations. 

Task 2. Construct logic for complete thermo-physical analysis 
within each physical device module. 

Task 3. Define physical device module interface with hardware 
performance characterization. 

Task 4. Construct specific computational logic sequencing 

routine for thermo-physical analysis of engine system. 

Task 5. Construct reconciler interface with device modules. 

Construct reconciler interface with gains model. 

Task 6. Write final report documenting integrated software 
platform development. 

SSME steady-state performance is defined by the fluid, flow 
and hardware characteristics that exist throughout the engine 
during steady-state operation. A logical platform for effective 
engine performance prediction must integrate physical principles 
and empirical information within a rigorous mathematical 
formalism. Physical principles pertinent to fluid flow and 
engine performance prediction are listed below. 
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Physical Principles 


1. Conservation of Mass 

2 . Conservation of Energy 

3 . The Second Law of Thermodynamics 

4 . Newton • s Second Law 

5 . Constitutive relations for thermal transport 
Hardware characteristics of particular interest in liquid-fueled 
rocket engine performance prediction include the following: 
Hardware Characteristics 

1. Turbine performance relations 

2. Pump performance relations 

3. Duct and other device flow resistance relations 

4. Chamber combustion efficiencies 

5. Nozzle efficiency. 

The objective of an efficient engine performance model is to 
identify physically consistent values of a complete, yet minimal, 
set of independent variables that describe the engine operating 
state. The variables selected must be consistent with the level 
of model approximation, and provide an adequate basis for 
insuring compatible hardware operating characteristics (see, e.g. 
[1] or [2]). 

For a one-dimensional steady-state flow model, the engine 
network can be defined by specifying the following information: 
Engine Network Components 

1. Flow branches, each consisting of a single inlet, 

single outlet flow path with associated flow rate and 
resistance 
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Engine Network Components (continued) 

2. Nodes, each defined as a specific location with an 
associated temperature and pressure 

3. Devices, each associated with a specific hardware 
component and defined by intersecting branches and 
boundary nodes 

4. Connectivity information defining each branch-branch 
and branch-device intersection and assigning individual 
nodes to each branch-branch and branch-device 
intersection . 

Notably absent from network definition components is detailed 
geometric information. Therefore, Newton's second law cannot be 
applied in a typical performance analysis to obtain fluid- 
structure interaction at the component level. 

As part of this investigation, a new one-dimensional steady- 
state rocket engine performance model has been constructed. This 
new model incorporates the physical principles described above, 
with the exception of Newton's second law, in a FORTRAN based 
computer software package. A simple interface for manufacturer 
supplied routines describing component hardware performance 
characteristics is provided. A description of the procedures 
employed in the new model is presented in the next section of 
this report. 
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2 . 0 ANALYSIS PROCEDURE 


The purpose of the new theoretical model is to provide an 
efficient computational procedure for defining engine flow 
networks and for determining physically consistent performance 
characteristics within the engine network. For SSME performance 
analysis, a minimal set of solution variables is prescribed 
below. Other operating characteristics can be derived from these 
primary variables. 

Primary SSME Performance Variables 

1. Mass flow rates through system branches [m (Ibm/s) ] 

2. Pressures at system nodes [P (psia) ] 

3. Temperatures at system nodes [T (R) ] 

4. Controllable valve resistances [R (s 2 /in 2 -ft 3 ) ] 

5. Turbopump shaft speeds [N (rpm) ] 

6. Heat transfer rates [Q (Btu/s) ] 

A number of computational strategies can be employed for 
one-dimensional steady-state network flow analysis. Adequate 
boundary and control setting specification is required to achieve 
solution closure. For SSME analysis, traditional inputs include 
flow inlet conditions (fluid composition, pressure, and 
temperature) , the desired oxygen-hydrogen mixture ratio, and the 
required thrust level or corresponding nozzle stagnation 
pressure. The following is a list of independent relations used 
to solve for the set of primary performance variables in the new 
engine model: 
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Model Analysis Relations 


1. Mass conservation at branch flow junctions 

2. Energy conservation in specified system volumes 

3. Pressure drop in flow branches as a function of branch 
resistance 

4 . Pressure drop across turbines as a function of turbine 
performance parameters 

5 . Temperature drop across turbines as a function of 
turbine performance parameters 

6. Head gain across pumps as a function of pump 
performance parameters 

7 . Power required by pumps as a function of pump 
performance parameters 

8 . System mixture ratio as a function of main combustion 
chamber (MCC) inlet flows 

9. Nozzle mass flow as a function of nozzle stagnation 
properties and geometry 

10. Heat transfer rate as a function of driving temperature 
difference and thermal resistance. 

These relations are presented below in residual form. An exact 

solution is a set of primary variable values that reduces each 

equation residual (or balance error) to zero: 

Governing Equations in Residual Form 

(mass flow residual); ( 1 ) 

i=l, 2 ,.. .number of I/O's in mass flow circuit j 
j=mass flow circuit number 



L my hjj - L m kj h kj + Qj = (energy flow residual )j (2) 

i k 

i=l , 2 ,... number of inputs to energy circuit j 
k=l, 2 ,.. .number of outlets from energy circuit j 
j=energy flow circuit number 
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Governing Equations in Residual Form ( cont inued ) 


( p m “ Pout “ R * W 2 / p)j = (pressure drop residual )j 
j=l, 2 , .number of pressure circuits 

{ p * N 2 * D 2 * [ C H - C H (C Q ) ] }j = (pump AP residual ) i 

j=l , 2 , . . . number of pumps 


[ 


[ 


m* ( P^-Pout) / P »* ( Pin-Pout) / P 


V 


3 = 1 , 2 ,. 


(Pin-Pout) -(Pin-? 


out/ characteristic- 

j ~ 1 r 2 f • 
(Pin - Pout) 


(Tj n ~T out ) “ (T^-T 


out/ characteristic - 

j = 1 r 2 , • 

(Pin-Pout) 


^ ( C Q / Ma ) J j 

. . number of pumps 


(pump power 
residual) j 


j = (turbine AP residual )j 
. . number of turbines 

characteristic Pin * ^ f / "Y ) 

j = (turbine AT residual )j 
. . number of turbines 

characteristic Pin * ^ ( ^'Q / / Y ) 


m oxygen/ m fuel " (^oxygen/^) command = (mixture ratio residual) 

E (nij) - m nozzJe = (nozzle flow residual) 
i 

i=l , 2 , ... number of nozzle inlet flows 
*W = f(P ns ,T ch , nozzle geometry) 


( Q - AT driving /R thcnnjd )j = (heat transfer residual^- 

j=l, 2 number of unknown heat transfer rates 

where 

m=mass flow rate 
h=specific enthalpy 
P=total pressure 
Q=heat transfer rate 
R=flow resistance 


(3) 

(4) 

(5) 

( 6 ) 

(7) 

( 8 ) 
(9) 

( 10 ) 
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T=temperature 

W=weight flow rate 

C Q =flow coefficient 

C H =head coefficient 

Ma=turbine Mach number 

P ns =nozzle stagnation pressure 

T ch =main combustion chamber temperature 

y=specific heat ratio 

i 7 =turbine efficiency 

p=mass density 

Inappropriate performance function 
These highly nonlinear relations are depicted graphically in 
Appendix B, Figures B2 through B6b. 

Many methods for solving systems of nonlinear equations are 
available (see, e.g. [3]). The objective of any solution 
procedure is to reduce the sum of all the residuals described 
above to zero. In practice this is generally not possible 
because of the approximate nature of both the physical relations 
and hardware performance curves. An appreciation of this 
limitation suggests the use of a minimization method to 
systematically and continuously reduce the residuals sum to an 
acceptable value. 

It is well known that the problem of solving a system of 
nonlinear equations may be replaced by a problem of minimizing a 
nonlinear function on R n [3], where n is the number of 
independent variables. In addition, the qua si -Newt on, or 
variable metric (VM) , methods are particularly robust strategies 
for minimizing unconstrained nonlinear functions. These facts 
motivated the use of a VM solution method in the new theoretical 
model. Specifically, a BFGS [4] multivariate search algorithm 
was implemented with an Armijo [5] univariate sub-algorithm. In 
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the present application, the objective of the BFGS-Armijo 
solution strategy was to select branch mass flow rates, nodal 
pressures and temperatures, variable valve resistances and 
turbopump shaft speeds in order to minimize the sum of the 
squares of the residuals defined above. To provide a consistent 
scale for the various residuals, each was divided by a user 
specified uncertainty estimate prior to squaring and summation. 
The engine operating point was thus obtained as the solution to 
the following mathematical programming problem: 


Minimize F - E 

i 


residual i 


uncertainty of residual i 
by selection of the primary performance variables. 


( 11 ) 


A stepwise description of the model analysis procedure is 
given below: 

Analysis Procedure 

1. Define network branches, nodes, and connectivity. 

2. Define mass, energy, pressure, turbopump, and nozzle 
circuits. 

3. Enter pertinent fluid property data. 

4 . Specify uncertainties including 

mass, energy, and pressure circuit balance 
uncertainties , 

turbine pressure and temperature drop 
uncertainties , 

pump pressure rise and power variance 
uncertainties , 
mixture ratio uncertainty, 
thrust uncertainty. 

5. Select desired mixture ratio and nozzle thrust. 
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Analysis Procedure (continued) 

6. Append turbopump performance curves. 

7. Append nozzle performance curves. 

8. Initialize branch flow rates, nodal pressures and 
temperatures, branch resistances, and turbopump speeds 
from available data. 

9. Solve for branch flow rates, nodal pressures and 
temperatures, valve controllable resistances, and 
turbopump speeds that provide engine balance. 

A hierarchy diagram displaying routines used in the new 
performance prediction model is presented in Figure Bl. A 
functional description of these routines together with a computer 
code listing of the new performance model is presented in 
Appendix C. Results based on preliminary testing of the new 
model are presented in the next section of the report. 
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3.0 PRELIMINARY RESULTS 


In order to test execution of the new model, a series of 
performance analyses was conducted on the SSME high pressure fuel 
turbopump subsystem depicted in Figure B7 . In this figure, nodes 
1 and 2 correspond to oxygen and hydrogen preburner inputs 
respectively. Nodes 3 and 4 are hot gas flow locations, and 
nodes 5 and 6 correspond to liquid hydrogen flow locations. Pump 
and turbine performance curves were approximated from predictions 
returned by Rocketdyne's SSME power balance model (PBM) over a 
typical range of engine power levels. All analyses were 
initiated at PBM solution values corresponding to 109% of SSME 
rated power level. Inlet temperatures and pressures at nodes 1, 
2, and 5 were fixed. In order to provide the basis for 
theoretical closure, command values of the preburner mixture 
ratio and pump power input were also designated at PBM solution 
values. 

Numerous test case analyses were conducted using somewhat 
arbitrary values of the residual scaling uncertainties. Results 
of three such analyses are presented in Appendix A. Table A1 
presents the scaling uncertainties used in each analysis. 

Specific values of these uncertainties were originally selected 
(Case A) to obtain common order of magnitude scaling of the 
various residual terms in Equation 11. In effect this made 
satisfaction of each residual relation approximately equal in 
importance . 
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In subsequent analyses, uncertainty scaling was used to 
emphasize or de-emphasize various relations. Case B uncertainty 
values were chosen identical to Case A values except for the 
preburner mass balance uncertainty which was an order of 
magnitude smaller in Case B. This indicates an increased 
emphasis on obtaining mass flow balance in Case B. Energy 
related uncertainties in Case C were reduced while uncertainties 
in the other turbine and pump performance relations were 
increased. This indicates an increased emphasis on obtaining 
overall power balance, and decreased confidence in hardware 
performance relations. 

Results of the three case analyses are presented in Tables 
A2 and A3 . Approximate solution values of the independent 
variables for each case are presented in Table A2 . Residuals 
associated with each of the governing balance relations are 
presented in Table A3 . In no case were all residual values 
reduced to zero identically. Therefore, an exact solution was 
not obtained in any of the three analyses. This was 
disconcerting although not surprising based on the approximate 
nature of the hardware performance curves used in the test cases, 
and the fact that the current PBM does not achieve exact 
solutions even with the best available performance information. 

Significant reductions in residuals were achieved at the 
approximate solutions returned by the new model. As shown in 
Table A3, the Case A solution provided significant reductions in 
power and pressure residuals accompanied by a small rise in 
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preburner flow residual. In order to reduce the flow residual 
(i.e., improve preburner mass flow balance), the flow balance 
uncertainty was reduced by an order of magnitude in Case B. A 
substantially different solution was obtained, with a decrease in 
flow residual from Case A and overall improvements in all power 
balances, although not as dramatic as provided by the Case A 
solution values. Similar comments could be made regarding the 
Case C solution residuals. 

The results of Case B were considered more realistic because 
of the low preburner flow imbalance. The Case B solution 
prescribes lower preburner and turbine flows, little change in 
initial pressure estimates, reduced turbine discharge 
temperature, and increased pump discharge temperature reflecting 
a decrease in pump efficiency. 
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4 . 0 RECOMMENDATIONS 


A list of recommendations based on construction, 
implementation, and run experience with the new performance model 
is presented below: 

1. Develop a strategy for assigning specific residual 
scaling uncertainty estimates for future model testing. 

2 . Interface existing turbopump performance curves to new 
model. 

3 . Expand current property routine to include potential 
engine states. 

4 . Streamline and structure property input interface to 
the new model. 

5. Implement postprocessing capability to recover 
additional hardware performance and design 
characteristics . 

6. Construct definition and connection input file for the 
full SSME engine network. 

7 . Perform an extensive computational test program on the 
new model applied to the full SSME engine system 
network in order to determine model efficiency and 
performance accuracy. 

8 . Compare new model computational results with existing 
power balance model predictions and Technology Test Bed 
(TTB) experimental data in order to assess integrity of 
existing PBM. 
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APPENDIX A 


TABLES 



Table Al. Summary of uncertainties for high pressure fuel turbopump 
subsystem analyses 


UNCERTAINTY ESTIMATE 
Preburner power balance (Btu/s) 
Preburner flow (lbm/s) 

0 2 resistance pressure drop (psi) 
H 2 resistance pressure drop (psi) 
Pump flow head gain (psi) 

Pump power requirement (hp) 
Turbine flow head drop (psi) 
Turbine flow temp drop (deg R) 
Turbopump power (Btu/s) 

Preburner 0 2 /H 2 ratio 


CASE A 

CASE B 

CASE C 

100 

100 

50 

H 

• 

O 

0.01 

0.01 

10 

10 

10 

10 

10 

10 

10 

10 

50 

100 

100 

50 

10 

10 

50 

1 

1 

10 

100 

100 

50 


0.001 


0.001 


0.001 



Table A2 


Summary of independent variable predictions for high 
pressure fuel turbopump subsystem analyses 


VARIABLE 

INITIAL 

CASE A 
FINAL 

CASE B 
FINAL 

CASE C 
FINAL 

Ml (lbm/s) 

79.63 

67.99 

74.68 

68.81 

M2 

84.08 

72.27 

78.12 

72.69 

M3 

163.72 

140.51 

152.80 

141.51 

M4 

162.25 

148.53 

152.16 

162.70 

PI (psia) 

7733.1 

(same) 

(same) 

(same) 

P2 

6088.0 

(same) 

(same) 

(same) 

P3 

5516.6 

5790.5 

5515.2 

5847.9 

P4 

3717.7 

3931.5 

3717.7 

3873.4 

P5 

271.7 

(same) 

(same) 

(same) 

P6 

6738.7 

6405.5 

6725.6 

6096.9 

T1 (deg R) 

208.3 

(same) 

(same) 

(same) 

T2 

278.4 

(same) 

(same) 

(same) 

T3 

1929.3 

1897.8 

1927.8 

1903.6 

T4 

1777.6 

1726.8 

1765.1 

1728.4 

T5 

42.4 

(same) 

(same) 

(same) 

T6 

96.7 

110.0 

105.0 

104.6 

N (rpm) 

36116 

34991 

36120 

36993 


Table A3 


Summary of imbalance predictions for high pressure fuel 
turbopump subsystem analyses 


IMBALANCE 

INITIAL 

CASE A 
FINAL 

CASE B 
FINAL 

CASE C 
FINAL 

Prebnr Flow 
(lbm/s) 

0.0004 

-0.2444 

0.002 

-0.0119 

Prebnr Power 
(Btu/s) 

-6493.7 

113.1 

-1265.2 

1009.1 

0 2 Res is AP 
(psia) 

-455.5 

-5.3 

-131.9 

-110.0 

H 2 Res is AP 
(psia) 

179.7 

8.1 

234.7 

-52.6 

HPFP AP 
(psia) 

-34.8 

-14.7 

-104 . 0 

-1015.6 

HPFP Power 
(hp) 

-10997.2 

-218.3 

-5209.2 

-1662.6 

HPFT AP 
(psia) 

12.7 

11.0 

8.9 

47.4 

HPFT AT 
(deg R) 

-0.0016 

24.6 

10.8 

22.1 

HPFTP Power 
(Btu/s) 

1006.1 

-794.7 

420.5 

-347.4 

0 2 /H 2 Ratio 

0.0 

-0.0063 

0.0089 

-0.0004 



APPENDIX B 


FIGURES 



FIGURE B1 


MODEL HIERARCHY DIAGRAM 
























FIGURE B3 • BRANCH PRESSURE LOSS 


R 



(P1-P2 ) - R * W**2 / Dens (PI, Tl) = PRESSURE RESIDUAL 


W = M * (gstd/gc) 



FIGURE B4 . BRANCH ENERGY BALANCE 


R 



M 


h(Pl,Tl) - h(P2,T2) = ENERGY RESIDUAL 


PI , P2 = TOTAL PRESSURES 



FIGURE B5a. 


TURBOPUMP BALANCES 


PI T1 P3 T3 



P2 T2 P4 T4 


Ml * [h (PI , Tl) - h (P2 , T2 ) ] + M2*[h(P3,T3) - h(P4,T4)] = POWER RESIDUAL 



FIGURE B5b. TURBOPUMP BALANCES 


Dens (PI , Tl) * (n*D) **2 * [CH - CH(CQ)] = PUMP DELTA P RESIDUAL 


CH 

CH(CQ) 

Eff 

Eff (CQ) 


= computed head coefficient 
= pump head coef as a function of flow coef 
= computed pump efficiency 

= pump efficiency as a function of flow coef 


Ml* (P2 - PI) /Dens (PI , Tl) 


Ml* (P2-P1) /Dens (PI , Tl) 


Eff 


Eff (CQ) 


PUMP POWER RESIDUAL 



FIGURE B5c 


TURBINE BALANCES 


[P3-P4 (assigned) -[P3-P4 (characteristic)] 
[T3-T4 (assigned) ] -[T3-T4 (characteristic)] 


TURBINE DELTA P RESIDUAL 
TURBINE DELTA T RESIDUAL 


P3-P4 

P3-P4 

(assigned) 

(characteristic) 

= program assigned delta P 

= delta P based on turbine characteristics 
= P3 * f [ CQ,Ma, gamma (P3,T3) ] 

T3-T4 

T3-T4 

(assigned) 

(characteristic) 

= program assigned delta T 

= delta T based on turbine characteristics 
= T3 * f [CQ,Ma, gamma (P3 ,T3) ] 


CQ 

Ma 

gamma (P3,T3) 


= turbine flow coefficient 
= turbine Mach number 

= hot gas specific heat ratio at P3 T3 



FIGURE B6a. NOZZLE AND OVERALL BALANCES 


PI T1 


P2 T2 



P3 T3 M3 


FIGURE B6b. NOZZLE AND OVERALL BALANCES 


M2/M1 - COMMANDED RATIO = MIXTURE RATIO RESIDUAL 


Ml + M2 - M3 = FLOW RESIDUAL 1 


M3 - M3 (P3,T3, nozzle geometry) = FLOW RESIDUAL 2 


Ml*h (PI , Tl) + M2*h (P2 , T2 ) - M3*h(P3,T3) = ENERGY RESIDUAL 


(P1-P3 ) - Rl*Wl**2/Dens(Pl / Tl) = PRESSURE RESIDUAL 1 


(P2-P3) - R2*W2**2/Dens (P2,T2) = PRESSURE RESIDUAL 2 


Thrust (M3 , P3 , T3 , nozzle geometry) - COMMAND THRUST = THRUST RESIDUAL 





APPENDIX C 


DOCUMENTATION 



Cl. DESCRIPTION OF ROUTINES 


ROUTINE 

MAIN 

START 

PROP 

PRPSAT 

PRPMIX 

ITERP2 

ITERP1 

OPT 

BFGS 

DFP 

ARMIJO 

OBJF 

GRAD 

TPIMB 

CURVES 


DESCRIPTION 

variable blocks configured, I/O data files identified, 
input data read 

solution variables initialized for solver, solver 
output routed to storage files 

property data arrays initialized, routing based on 
input property values provided 

NBS properties of parahydrogen and oxygen near the 
respective saturation curves calculated 

properties of gaseous mixtures of parahydrogen, oxygen, 
and water calculated using the Dalton model 

table constructed functions of two variables 
interpolated to estimate function value at non-table 
independent variable values 

table constructed functions of one variable 
interpolated to estimate function value at non-table 
independent variable values 

iterative logic sequence for BFGS-Armijo solver 
provided 

BFGS variable metric Hessian update calculated 

Davidon-Fletcher-Powell variable metric Hessian update 
calculated 

Armijo univariate search conducted to determine 
residuals minimum in search direction provided by OPT 

residuals function of system governing equations 
calculated 

forward difference finite difference approximation of 
residuals function gradient calculated 

turbine and pump residuals based on turbopump 
performance curves calculated 

turbine and pump performance characteristics provided 
by manufacturer 



C2 • INPUT/OUTPUT VARIABLE DEFINITIONS 


Program Input File Name = TM2-I0.DAT 
Data Type = I/O data file identification 
INPUT VARIABLE NAME DEFINITION 


PDAT 

TDAT 

UDAT 

VDAT 


name of file containing fluid property data 
name of file containing solution initiation data 
name of file containing residual uncertainty 
estimates 

name of file containing volume definition and 
connection information 

name of file containing solution output 


ODAT 



Program Input File Name 


VDAT 


Actual Access File Name = "named in TM2-I0.DAT” 

Data Type = system component definition and connection data 


INPUT VARIABLE NAME 


DEFINITION 


NEC 

NMC 

NPC 

NTPC 

NBRH 

NNOD 

NHGNOD 

NRPM 

INOZBR 

INOZEC 

INOZMC 

INOZN 

IHGNOZ 

IMIXR 

IATHRT 

IA (I) 
IM(I) 

IP(I) 

IR(I) 

IS(I) 

IT (I ) 

MAT (I) 


NEIO(I) 
EDIR( I , J) - 


IEN (I , J) - 
IMBEC (I , J) - 
NMIO(I) 
MDIR(I,J) - 


IMBMC ( I , J ) - 


number of energy circuits in system 

number of mass flow circuits in system 

number of pressure circuits in system 

number of turbopump circuits in system 

number of mass flow branches in system 

number of energy (P,T) nodes in system 

number of hot gas nodes in system 

number of turbopump shaft speeds 

mass flow branch number of nozzle flow 

energy circuit number containing nozzle flow 

mass flow circuit number containing nozzle flow 

energy node number of nozzle flow 

hot gas node number of nozzle flow 

TTB array location containing mixture ratio value 

TTB array location containing nozzle throat area 

value 

(not used) 

TTB array location containing initial estimate of 
mass flow rate through branch I 

TTB array location containing initial estimate of 
pressure value at node I 

TTB array location containing initial estimate of 
flow resistance value for branch I 
TTB array location containing initial estimate of 
turbopump I shaft speed 

TTB array location containing initial estimate of 

temperature value at node I 

material type number of flow at node I 

1 = hydrogen (H2) 2 = oxygen (02) 

3+= hot gas mixture 

number of mass flow I/O's across boundary of 
energy circuit I 

direction number of mass flow J across boundary of 
energy circuit I 

1 = inflow -1 = outflow 

node number associated with flow J across boundary 
of energy circuit I 

branch number of mass flow J across boundary of 
energy circuit I 

number of mass flow I/O's across boundary of mass 
circuit I 

direction number of mass flow J across boundary of 
mass circuit I 

1 = inflow -1 = outflow 

branch number of mass flow J across boundary of 
mass circuit I 



VDAT (continued 2) 


Program Input File Name = 

INPUT VARIABLE NAME DEFINITION 

PDIR(I,J) - direction number of mass flow J across boundary of 

pressure circuit I 

1 = inflow -I = outflow 

IPN(I,J) - node number of mass flow J across boundary of 

pressure circuit I 

IMBPC(I,J)- branch number of mass flow J across boundary of 

pressure circuit I 

NH2HG (I) - number of pure hydrogen flows contributing to hot 

gas flow I 

N02HG(I) - number of pure oxygen flows contributing to hot 

gas flow I 

IHGN(I) - node number associated with hot gas flow I 
ICEFF(I) - TTB array location containing the combustion 

efficiency value associated with the flow at hot 
gas node I 

IH2HG(I,J)- branch number of hydrogen flow J entering hot gas 

flow I 

I02HG(I,J)- branch number of oxygen flow J entering hot gas 

flow I 

IBTYPE(I) - type number of branch flow I 

0 = fixed >0 = variable 

INTYPE (I) - type number of node I 

0 = fixed pressure and temperature 

1 - variable pressure, fixed temperature 
>1 = variable pressure and temperature 

IRTYPE (I) - type number of resistance in pressure circuit I 

0 = fixed resistance 

>0 = variable resistance 
ITPTY(I) - type number of turbopump I 

1 = one turbine, one pump 
>1 = one turbine, two pumps 

ITPD (I , J) - TTB array location containing the impeller 

diameter of turbomachine J in turbopump circuit I 
J = 1 turbine J = 2 or 3 pump 

ITPN(I,J) - node number associated with mass flow J crossing 

boundary of turbopump I 

J = 1 turbine inlet J = 2 turbine outlet 

J = 3 pump 1 inlet J = 4 pump 1 outlet 

J = 5 pump 2 inlet J = 6 pump 2 outlet 

ITPS(I) - shaft number associated with turbopump I 
IMBTP (I , J) - branch number associated with mass flow J crossing 

boundary of turbopump I 

J = 1 turbine inlet 

J = 3 pump 1 inlet 

J = 5 pump 2 inlet 


J = 2 turbine outlet 
J = 4 pump 1 outlet 
J = 6 pump 2 outlet 



VDAT (continued 3) 


Program Input File Name = 

INPUT VARIABLE NAME DEFINITION 

ICURV(I,J)- curve number associated with performance curve J 

of turbopump I 

J = 1 turbine performance curve 1 
J = 2 turbine performance curve 2 
J — 3 pump 1 performance curve 1 

J — 4 pump 1 performance curve 2 

J = 5 pump 2 performance curve 1 

J = 6 pump 2 performance curve 2 



Program Input File Name = TDAT 

Actual Access File Name = "named in TM2-I0.DAT" 

Data Type = analysis description and solution search 
initializing values 


INPUT VARIABLE NAME DEFINITION 


NDESC 

NTTB 

DESC(I) 

TTB(I) 


number of analysis description items 
number of elements in the TTB solution 
initializing array 

analysis description item I (character variable) 
value of solution initializing variable assigned 
to array address I 


Program Input File Name = UDAT 

Actual Access File Name = "named in TM2-IO.DAT" 

Data Type = uncertainty values associated with model delations 

INPUT VARIABLE NAME DEFINITION 


UOF 

UEVOL(I) 

UMVOL(I) 

UPVOL(I) 

UETP(I) 

UTITP(I) 

UT2TP (I) 

UPITP(I) 
UP2TP (I) 
UP3TP (I) 
UP4TP (I) 


uncertainty of engine oxygen/fuel ratio 

power imbalance uncertainty associated with energy 

circuit I 

mass flow rate imbalance uncertainty associated 
with mass flow circuit I 

pressure imbalance uncertainty associated with 
pressure circuit I 

power imbalance uncertainty associated with 
turbopump circuit I 

turbine characteristic 1 (pressure drop) imbalance 
uncertainty associated with turbopump circuit I 
turbine characteristic 2 (temperature drop) 
imbalance uncertainty associated with turbopump 
circuit I 

pump 1 characteristic 1 (pressure drop) imbalance 
uncertainty associated with turbopump circuit I 
pump 1 characteristic 2 (power) imbalance 
uncertainty associated with turbopump circuit I 
pump 2 characteristic 1 (pressure drop) imbalance 
uncertainty associated with turbopump circuit I 
pump 2 characteristic 2 (power) imbalance 
uncertainty associated with turbopump circuit I 



Program Output File Name 


ODAT 


Actual Access File Name = "named in TM2-IO.DAT" 

Data Type = output defining circuit definitions, solution 
values, and residuals (both initial and final) 


INPUT VARIABLE NAME 


DEFINITION 


NEC 
NMC 
NPC 
NTPC 
IEN (I , J) 

IMBEC(I, J)- 

IMBMC ( I , J ) - 

IPN (I , J) - 

IMBPC (I , J) - 

ITPN (I , J) - 

IMBTP (I , J) - 

ETPO(I) 
ETP(I) 
TITPO(I) - 
T1TP(I) 
T2TP0 (I) - 

T2TP (I) 
PITPO(I) - 
P1TP(I) 
P2TP0 (I) - 

P2TP (I) 


number of energy circuits in system 

number of mass flow circuits in system 

number of pressure circuits in system 

number of turbopump circuits in system 

node number associated with flow J across boundary 

of energy circuit I 

branch number of mass flow J across boundary of 
energy circuit I 

branch number of mass flow J across boundary of 
mass circuit I 

node number of mass flow J across boundary of 
pressure circuit I 

branch number of mass flow J across boundary of 
pressure circuit I 

node number associated with mass flow J crossing 
boundary of turbopump I 

J = l turbine inlet J = 2 turbine outlet 
J = 3 pump 1 inlet J = 4 pump 1 outlet 

J = 5 pump 2 inlet J = 6 pump 2 outlet 

branch number associated with mass flow J crossing 
boundary of turbopump I 

J = 1 turbine inlet J = 2 turbine outlet 

J = 3 pump 1 inlet J - 4 pump 1 outlet 

J = 5 pump 2 inlet J = 6 pump 2 outlet 

initial power imbalance associated with turbopump 
circuit I 

final power imbalance associated with turbopump 
circuit I 

initial turbine characteristic 1 (pressure drop) 
imbalance associated with turbopump circuit I 
final turbine characteristic 1 (pressure drop) 
imbalance associated with turbopump circuit I 
initial turbine characteristic 2 (temperature 
drop) imbalance associated with turbopump circuit 
I 

final turbine characteristic 2 (temperature drop) 
imbalance associated with turbopump circuit I 
initial pump 1 characteristic 1 (pressure drop) 
imbalance associated with turbopump circuit I 
final pump 1 characteristic 1 (pressure drop) 
imbalance associated with turbopump circuit I 
initial pump 1 characteristic 2 (power) imbalance 
associated with turbopump circuit I 
final pump 1 characteristic 2 (power) imbalance 
associated with turbopump circuit I 



Program Output File Name 
INPUT VARIABLE NAME 


ODAT (continued 2) 
DEFINITION 


P3TP0 (I) 
P3TP (I) 
P4TP0 (I) 
P4TP (I) 


initial pump 2 characteristic 1 (pressure drop) 
imbalance associated with turbopump circuit I 
final pump 2 characteristic 1 (pressure drop) 
imbalance associated with turbopump circuit I 
initial pump 2 characteristic 2 (power) imbalance 
associated with turbopump circuit I 
final pump 2 characteristic 2 (power) imbalance 
associated with turbopump circuit I 



C3 • SOURCE CODE LISTING 





PROGRAM TM-02 


C 

C 


c 


CHARACTER* 2 4 DESC 

CHARACTER* 12 PDAT , TDAT , UDAT , VDAT , ODAT 
INTEGER EDIR, PDIR 


COMMON /BALC/ EIMB(5) 
EIMBO (5) 
ETP (4 ) 
P1TP (4 ) 
ETPO (4 ) 
P1TP0 (4 ) 
RE VA (20) 
REVR (20) 


, PIMB ( 6 ) 

, PIMBO ( 6 ) 
, T1TP (4 ) 

, P2TP ( 4 ) 
TITPO (4 ) 
, P2TP0 (4 ) 
, RE VM (20) 
, REVS (20) 


WIMB (5) 
WIMBO (5) 

, WZIMB 

OFIMB 

, WZIMBO , 

OFIMBO 

T2TP I 
P3TPI 

[t] 

J P4TP ( 4 ) , 


T2TPC 


/ 


P3TP0 (4 ) 

, P4TP0(4), 


REVPI 

[ 2 ° 

, REVT(20) , 


REVDI 

[ 20 ) 

, REVH (20) 



C 

c 


c 


c 


c 


c 

c 


c 


c 


COMMON /VDAT/ NEC, NMC, NPC, NTPC, NBRH, NNOD, NHGNOD , NRPM, 
INOZBR, INOZEC , INOZMC , INOZN , IHGNOZ , IMIXR, I ATHRT , 
IA(20) . IM(20) , IP (20) ,IR( 6 ) ,IS(5) , IT (20) ,MAT (20) , 


NElO(5 
NMIO (5) 
IPN($,2) 
IHGN (5) 
IBTYPE (20) 
ITPD (4,3) 
ICURV (4,6) 


EDIR (5 ,20) , 
MDIR( 5 ,20) , 
IMBPC ( 6 ) , 

ICEFFX5) 
INTYPE (20) , 
ITPN (4,6) , 


±EN (5,20) 
IMBMC(5,20) 
NH2HG (5) 
IH2HG|5,5) 
IRTYPE(20) 
ITPS (4 ) 


IMBEC (5,20) 
PDIR (6.2) 
N02HG (5) 


I02HG 

ITPTY 

IMBTP 


Ij 5) 

4,3) 


COMMON /TDAT/ NDESC , NTTB, DESC(5), TTB(1350) 


COMMON /UDAT/ UOF , UEVOL(5) , UMVOL(5 

1 UETP ( 4 ) , UT1TP (4 ) , UT2TP(4' 

3 UP2TP (4 ) , UP3TP (4 ) , UP4TP(4 


UPVOL ( 6 ) , 
UP1TP ( 4 ) , 


COMMON /H2PRP/ 

* H2P1 (15) , H2T1 ( 11) , H2H1 (15,11) ,H2S1(15,11) ,H2D1(15,11 

* H2P2 (20) ,H2T2 (11) ,H2H2 (20,11 ,H2S2 (20,11 ,H2D2 (20,11 

* H2P3 (29) , H2T3 (25 ,H2H3 (29,25 ,H2S3 (29,25 ,H2D3 (29,25 

* H2P4 (23 ) , H2T4 (25) ,H2H4 (23,25 ,H2S4 (23,25) ,H2D4(23,25 
COMMON /02PRP/ 

* 02P1 ( 13 ) , 02T1 ( 16 ) , 02H1 (13,16) ,02S1(13,16) ,02D1(13,16) , 

* 02P2 ( 13 ) , 02T2 ( 17 j ,02H2(13,17) ,02S2 (13,17) ,02D2(13,17) , 

* 02P3 (5) , 02T3 (61) ,02H3 (5, 61 , 02S3(5,61), 02D3(5,61) 
COMMON / 62 OPRP/ 

* H20P1 (7) , H20T1 ( 13 ) ,H20H1(7,13) ,H20S1(7,13) ,H20D1(7,13) 
COMMON /TABLE/ 

* NH2P (4 ) ,NH2T (4 ) ,N02P(3) ,N02T(3) ,NH20P(1) ,NH20T(1) 
COMMON /STD/ 

* HH2REF , H02REF , HWAREF , SH2REF , S02REF , SWAREF , SH2A , S02 A , 

* SWA A 


DIMENSION 

* NH2PA ( 4 ) , NH2TA ( 4 ) ,N02PA(3) ,N02TA(3) ,NH20PA(1) ,NH20TA(1) 
CHARACTER* 70 PTITLE 


DATA (NH2PA (I ) , 1=1 , 4 ) /15 , 20 , 29 , 23/ 
DATA (NH2TA ( J) , J=1 , 4 ) /II , 11 , 25 , 25/ 
DATA (N02PA(I) , 1=1, 3)/13,13,5/ 

DATA (N02TA j[ , J=1 , 3 ) /16 , 17 , 61/ 
DATA NH20PA(I) ,1=1,1) /7/ 

DATA (NH20TA(J) ,J=1,1)/13/ 

DO 90 1=1,4 
NH2 P ( I ) — NH2 PA ( I ) 

90 NH2T (I J =NH2TA(I) 

DO 91 1=1,3 

N02 P ( I ) =N02 PA ( I ) 

91 N02T]i)=N02TA(I) 

NH20P (1) =NH20PA (1) 

NH20T(1 =NH20TA(1) 


HH2REF 

H02REF 

HWAREF 

SH2REF 

S02REF 

SWAREF 

SH2A 


1790.091 

234.681 

1339.990 

15.440 

1.530 

2.294 

15.481 
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up rso H 'O CT\ ttv U> H 


C 

c 


REAL JOULE 

DIMENSION X(25) ,X2 (25) 

COMMON /BALC/ EIMB(5) , PIMB(6) , WIMB(5) , WZIMB , OFIMB , 

EIMB0(5) , PIMBO (6) , WIMB0(5) , WZIMBO , OFIMBO, 

ETP (4 ) , T1TP(4) , T2TP (4 ) 

P1TP (4 ) , P2TP (4 ) , P3TP (4 ) 

ETPO (4) , T1TP0 (4 ) , T2TP0(4 

P1TP0 (4) , P2TP0 ( 4 ) , P3TP0 (4 

REVA ( 2 0 j , RE VM (20) , REVP(20 

REVR(20) , REVS (20) , REVD(20 


' P4TP (4 ) , 

P4TP0 (4) , 
, RE VT (20) , 
, REVH (20) 


COMMON /VDAT/ NEC, NMC, NPC, NTPC, NBRH, NNOD, NHGNOD , NRPM, 
INOZBR, INOZEC, INOZMC, INOZN, IHGNOZ , IMIXR, IATHRT 
IA(20) , IM(20) , IP(20) , IR(6) , IS (5) , IT (20) ,MAT (20 
NEIO(5) , EDIR(5, 20) , IEN(5,20) , IMBEC(5,: 

NMIO (5) . MDIR (5,20) , IMBMC(5,20), PDIR(6,2 


IA ( 2 
NEIO 
NMIO 
IPN ( 
IHGN 
IBTY 


lil! 


IBTYPE (20) 
ITPD C4 , 3 ) 
ICURV (4,6) 


ICEFF [ 5 I 
INTYPE (20 
ITPN (4,6) 


IH2HG (5 
IRTYPE ( 
ITPS (4) 


!sS) 


I02HG (5.5) 
ITPTY ( 4 ) 
IMBTP (4,3) 


/ 

/ 

/ 

/ 

/ 

/ 


COMMON /TDAT/ NDESC, NTTB, DESC(5), TTB(1350) 


COMMON /UDAT/ UOF , UEVOL(5) , UMVOL(5) , UPVOL(6), 

UETPj4) , UT1TP(4) , UT2TP (4 ) , UP1TP(4), 

UP2TP ( 4 ) , UP3TP(4) , UP4TP(4) 

C 

COMMON /H2PRP/ 

1 H2P1 (15) ,H2T1 ( 11) , H2H1 (15,11) ,H2S1(15,11) ,H2D1(15,11) , 

2 H2P2 (20) ,H2T2 (11) ,H2H2 (20,11) ,H2S2 (20,11) ,H2D2 (20,11), 

3 H2P3 (29) , H2T3 (25) ,H2H3 (29,25) ,H2S3 (29,25) ,H2D3 (29,25) , 

4 H2P4(23 ,H2T4 (25) ,H2H4 (23,25) ,H2S4(23,25) ,H2D4 (23,25) 

COMMON /02PRP/ 

1 02P1 ( 13 ) , 02T1 ( 16 ) , 02H1 (13,16) ,02S1(13,16) ,02D1(13,16) , 

2 02P2 ( 13 ) , 02T2 ( 17 ) ,02H2 (13,17) ,02S2(13,17) ,02D2(13,17) , 

3 02P3 (5) , 02T3 (61) ,02H3 (5, 61) , 02S3(5,61), 02D3(5,61) 

COMMON /H20PRP/ 

1 H20P1 (7) ,H20T1 (13 ) ,H20H1(7,13) ,H20S1(7,13) ,H20D1(7,13) 

C 

COMMON /STD/ 

1 HH2REF , H02REF , HWAREF , SH2REF , S02REF , SWAREF , SH2A , S02 A , 

2 SWAA 

COMMON /TABLE/ 

1 NH2P(4) , NH2T ( 4 ) ,N02P(3) ,N02T(3) ,NH20P(1) ,NH20T(1) 

C 

PARAMETER ( JOULE = 778.16, GC = 32.174 ) 

C 

IFUNC=0 

C 

WRITE ( 21, 941 ) 

DO 2 IMC = 1, NMC 

WRITE ( 21, 942 ) IMC, ( IMBMC ( IMC , IIO) , IIO = 1, NMIO ( IMC ) ) 
2 CONTINUE 
C 

WRITE ( 21, 943 ) 

DO 4 IEC = 1, NEC 

WRITE ( 21, 942 ) IEC, ( IMBEC (IEC, IIO) , IIO = 1, NEIO ( IEC ) ) 

WRITE ( 21, 944 ) IEC, ( IEN (IEC, IIO) , IIO = 1, NEIO ( IEC ) ) 

4 CONTINUE 
C 

WRITE ( 21, 945 ) 

DO 6 IPC = 1, NPC 

WRITE ( 21, 942 ) IPC, IMBPC(IPC), IMBPC(IPC) 

WRITE ( 21, 944 ) IPC, ( IPN (IPC, IIO) , IIO = 1, 2 ) 

6 CONTINUE 
C 

WRITE ( 21, 946 ) 

DO 7 ITPC = 1, NTPC 


WRITE | 

[ 21, 

942 ) 

ITPC, 

( IMBTP 1 

[ITPC, J) , 

J = 

1, 

3 ) 

WRITE I 
7 CONTINt 

[ 21, 
JE 

944 ) 

ITPC, 

( ITPN I 

[ITPC, J) , 

J = 

1/ 

6 j 


C 

C ********** INITIALIZE ALGORITHM ************************************* 
C 


NMAX=100 



ooo 


EPSC=1 . E-6 
EPST=1 . E-6 
RH0=1 . 0 
ALPHA=Q . 0 
BETA=0 . 4 

********** INITIALIZE INDEPENDENT VARIABLES X(I) 


10 


c 


20 


c 


40 


c 


50 

c 

c 


c 


c 


c 


60 


70 

72 


K=1 

DO 10 I=1,NBRH 
RE VM? I 1=TTB ( IM ( I ) ) 

IF (IBTYPE(I) .EQ.O) THEN 
GO TO 10 
ELSE 

X^K) =SQRT (REVM ( I ) ) 

ENDIF 

CONTINUE 

DO 20 1=1 , NNOD 
REVP ( I ) =TTB ( IP ( I ) ) 

REVT? I) =TTB (IT (I) 

IF (INTYPE (I) .EQ.O) THEN 
GO TO 20 
ELSE 

IF (I.EQ.INOZN) THEN 
a (K) =SQRT (REVT (I) ) 
K=K+1 
ELSE 

X (K) =SQRT (REVP (I) ) 

X ( K4- 1 ) =SQRT ( REVT ( I ) ) 
K=K+2 
ENDIF 
ENDIF 
CONTINUE 

DO 40 1=1, NPC 
REVR ( I )=TTB ( IR (I ) ) 

IF (IRTYPE(I) .EQ.O) THEN 
GO TO 40 
ELSE 

X_(K) =SQRT (REVR (I) ) 

ENDIF 

CONTINUE 

DO 50 1=1 , NRPM 
REVS (I ) =TTB ( IS ( I ) ) 

X (K) =SQRT (REVS (I) ) 

K=K+1 

CONTINUE 

N=K-1 


DO 100 1=1, NNOD 

P=REVP( I ) 

T=REVT( I ) 


IF i 

f MAT (I 

) .GE. 

3 ) 

GO 

TO 

70 

IF i 

( MAT (I 

) .GE. 

2 ) 

GO 

TO 

60 


CALL PROP ( 1, P, T, 0.0, 0.0, H, S, RHO) 
REVD ( I ) =RHO 
REVH ( I ) =H-HH2REF 
GO TO 100 

CALL PROP ( 2, P, T, 0.0, 0.0, H, S, RHO) 
REVD ( I ) =RHO 
REVH(I) == H-H02REF 
GO TO 100 

IDHG=1 

IF (IHGN(IDHG) .EQ.I) THEN 
IHG=IDHG 
GO TO 74 
ELSE 

IDHG=IDHG+1 

IF (IDHG.GT.NHGNOD) THEN 


******************** 



c 


74 


80 


C 


90 


C 


c 

100 


GO TO 74 
ELSE 

GO TO 72 
ENDIF 
ENDIF 


WH2=0 . 0 

DO 80 IH2IN=1 , NH2HG (IHG) 
IBRH=IH2HG ( IHG , IH2 IN) 
WH2=WH2+REVM ( IBRH) 
CONTINUE 


W02=0.0 

DO 90 I02IN=1 , N02HG (IHG) 
IBRH=I02HG (IHG. I02IN) 

WO 2 =WO 2 +RE VM ( I BRH ) 
CONTINUE 


CEFF=TTB (ICEFF ( IHG ) ) 

0F=W02/WH2 

CALL PROP ( 4, P, T, OF, CEFF, H, 
REVD ( I ) =RHO 
REVH ( I ) =H 

CONTINUE 


C 

c ********** CALL OPTIMIZATION STRATEGY 
C 


S, RHO) 


******************************* 


CALL OPT (N,NMAX,EPSC.EPST, RHO, ALPHA, BETA, ISTAGE , F2 , GNORM, 
1 DXNORM , X , X2 ) 

C 

WRITE (21,801) 

DO 200 I=1,NBRH 

WRITE (21,802) I , TTB ( IM ( I ) ) , REVM ( I ) 

200 CONTINUE 
C 

DO 210 1=1 ,NNOD 

WRITE (21,803) I , TTB ( IP ( I ) ) , REVP ( I ) 

210 CONTINUE 
C 


220 

C 


230 

C 


240 

C 


250 

C 


260 

C 


270 


DO 220 1=1 , NNOD 

WRITE (21,804) I , TTB (IT (I) ) ,REVT (I) 
CONTINUE 


DO 230 1=1 , NRPM 

WRITE (21,805) I , TTB (IS ( I ) ) , REVS (I) 
CONTINUE 

WRITE (21,806) 

DO 240 1=1, NMC 

WRITE (21,807) I , WIMB0 (I) , WIMB (I) 
CONTINUE 

DO 250 1=1, NEC 

WRITE (21,808) I , EIMBO ( I ) , EIMB ( I ) 
CONTINUE 


DO 260 1=1, NPC 

WRITE (21,809) I, PIMBO (I) , PIMB(I) 
CONTINUE 


1 
2 

3 

4 

ELSE 

WRITE (21,811) 

1 

2 

3 

4 

5 

6 

ENDIF 

CONTINUE 


, ETP(I) 
, T1TP (I 
, T2TP (I 
, P1TP (I 
, P2TP ( I 

I,ETP0(I) , ETP (I) 
I , T1TP0 (I) ,T1TP(I 
I,T2TP0 (I) , T2TP ( I 
I , P1TP0 (I) , P1TP ( I 
I , P2TP0 (I) , P2TP ( I 
I , P3TP0 (I) , P3TP ( I 
I , P4TP0 (I) , P4TP ( I 


DO 270 1=1, NTPC 
IF (ITPTY ( I ) . EQ . 1 ) THEN 

WRITE (21,810) I , ETPO (I ) 

I , T1TP0 (I 
I , T2TP0 ( I 
I , P1TP0 ( I 
I , P2TP0 (I 


/ 


/ 



c 


WRITE (21,812) I , OFIMBO , OFIMB 
WRITE (21,813) I , WZ IMBO , WZ IMB 


801 FORMAT 


802 

803 

804 

805 

806 

807 

808 

809 

810 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


1 
2 

3 

4 

811 FORMAT 
1 
2 

3 

4 

5 

6 


(/12X 'ANALYSIS 
'INITIAL' , 15X, 
IX,' FLOW', 

IX, 'PRESSURE' , 
IX,' TEMP', 

(IX,' SPEED' 
/16X, 'BALANCES 
'FINAL' ) 

5X, ' MASS ' , 

5X, ' ENERGY', 
‘5X. 'PRESSURE' , 
; IX, ' ENERGY 
/IX, ' TURB DP 
/IX, ' TURB DT 
/IX, ' PUMP DP 
/IX, 'PUMP PWR 


RESULTS' ,/ 4X, 'VARIABLE' ,15X, 
' FINAL ' ) 

I3,7X,F15.4,5X, 

13 , 7X, FI 5 . 4 , 5X, 

13 , 7X, F15 . 4 , 5X, 


F15 . 4 , 5X, ' (LB/S) ' ) 
F15 . 4 , 5X, ' (PSI) * ) 
F15 . 4 , 5X, ' (DEG R) ' ) 
13 , 7X, F15 . 4 , 5X, F15 . 4 . 5X, ' XRPM) ' ) 

' ,/8X, ' CIRCUIT ' , 12X, ' INITIAL' ,12X, 


( ix, 

/IX, 

/IX, 

/IX, 

/IX, 

/IX, 

/IX, 


12 , 7X, F12 . 4 , 5X, 
12 , 7X, F12 . 4 , 5X, 
I2,7X,F12 .4 ,5X, 
STP ' , 12 , 7X , F12 . 
STP ' , 12 , 7X, F12 . 
STP' , 12 , 7X, F12 . 
STP' , 12 , 7X, F12 . 
STP' , 12 , 7X, F15 . 


F12 . 4 , 5X, ' (LB/S) ' ) 


ENERGY 
TURB DP 
TURB DT 
PUMP1 DP 
PUMP1 PW 
PUMP 2 DP 
PUMP 2 PW 


BTP 

BTP 

BTP 

BTP 

BTP 

BTP 

BTP 



4 , 5X, F12 
4 , 5X, F12 . 4 , 5X, 
4 , 5X, F12 . 4 , 5X, 
4 , 5X, F12 . 4 , 5X, 
4 , 5X, F15 . 4 , 5X, 


, 12 , 7X, F12 . 4 , 5X , F12 . 4 , 5X , 
, 12 , 7X, F12 . 4 , 5X, F12 . 4 , 5X, 
, 12 , 7X , F12 . 4 , 5X , F12 . 4 , 5X , 
, 12 , 7X, F12 . 4 , 5X, F12 . 4 , 5X, 
, 12 , 7X , F15 . 4 , 5X, F15 . 4 , 5X , 
, 12 , 7X, F12 . 4 , 5X, F12 . 4 , 5X, 
,I2,7X,F15.4, 5X, F15 . 4 , 5X, 


TU/S) ' , 
PSI) ' , 

r T 1 # 

PSI) ' . 
BTU/S) ' ) 


BTU/S) 

Rvl ' ’ 

Mu 


) 


812 

813 

FORMAT 

1 

FORMAT 

1 

(/ix, • 

F7.3) 

941 

FORMAT 

1 

( >/4X, 

942 

943 

FORMAT 

FORMAT 

( 8X, 

( / t 2X f 

944 

945 

FORMAT 

FORMAT 

( 8X , 

( /fix, 

946 

FORMAT 

1 

( /f9X, 


RETURN 

END 



O/F - O/F COMMANDED' ,5X, 'INITIAL' ,F7. 3, 5X, ’FINAL' , 
NOZZLE FLOW IMBALNC ', 5X ,' INITIAL ', F7 . 3 , 5X, ' FINAL ' , 


, 5X, 'I/O NUMBER', 6X, 
5 6' ) 


'FLOW SUBSYSTEM NO' 

12 3 4 

5, 12X, 'I/O FLOW BRANCHES 615 ) 
'ENERGY SUBSYSTEM NO', 5X, ' I/O NUMBER 

'■ ) 


..*.3 4 5 6 ' 

!5, 12X, 'I/O ENERGY NODES ','615 
'DELTA P SUBSYSTEM NO', 5X, 'I/O 
1 2 ‘ 




UMBER 


PU^] 


6X, 


6X, 


i 


********************************************************************* 


SUBROUTINE OBJF (IFUNC,X,F) 

CHARACTER* 2 4 DESC 
INTEGER EDIR, PDIR 

REAL JOULE 


DIMENSION X (25) 

COMMON /BALC/ EIMB(5) 

1 EIMBO (5) 

3 ETP(4) 

4 P1TP (4 

6 ETPO ( 4 

7 P1TP0 ( 4 

1 REVA (20 

2 REVR (20 


PIMB (6) 
PIMBO ( 6) 
T1TP (4 ) 
P2TP ( 4 ) 
T1TP0 (4 
P2TP0 ( 4 
REVM (20 
REVS (20 


WIMB ( 5 ) 
WIMBO (5) 
T2TP (4 ) 
P3TP ( 4 ) 
T2TP0 (4 
P3TP0 ( 4 
REVP (20' 
REVD ( 2 0 


WZIMB 

WZIMBO 

P4TP ( 4 ) 

P4TP0 ( 4 
REVT (20 
REVH (20 


OFIMB , 
OFIMBO, 


COMMON /VDAT/ NEC, NMC, NPC, NTPC, NBRH, NNOD, NHGNOD, NRPM, 
INOZBR, INOZEC, INOZMC, INOZN, IHGNOZ , IMIXR, IATHRT. 
IA (20) , IM (20) ,IP?20) , IR ( 6 ) , IS ( 5 ) , IT (20) , MAT (20) , 


NEIOl 
NMIO (5) 

IPN (6,2) 
IHGN (5) 
IBTYPE (20) 
ITPD (4,3) 

I CUR V (4,6) 


EDIR(5, 20 
MDIR (5,20 
IMBPC ( 6 ) 
ICEFF (5) 
INTYPE (20) 
ITPN (4,6) 


±EN (5,20) , IMBEC 
IMBMC (5,20) , PDIR ( 
NH2HG (5) , N02HG 
IH2HG? 5 , 5 ) , I02HG 
IRTYPE (20) , ITPTY 
ITPS ( 4 ) , IMBTP 




c 


c 


c 


c 


c 

c 

c 


10 


c 


20 


c 


40 


c 


50 


c 


c 


c 


COMMON /TDAT/ NDESC, NTTB, DESC(5), TTB(1350) 

COMMON /UDAT/ UOF , UEVOL(5) , UMVOL(5) , UPVOL(6) , 
1 UETP ( 4 ) , UT1TP (4 ) , UT2TP(4), UP1TP(4), 
3 UP2TP ( 4 ) , UP3TP ( 4 ) , UP4TP(4) 


COMMON /H2PRP/ 

1 H2P1(15) ,H2T1(11) ,H2H1(15,11 

2 H2P2 (20) ,H2T2 (11) ,H2H2 (20, 11 

3 H2P3 (29) ,H2T3 (25) ,H2H3 (29,25 

4 H2P4 (23) , H2T4 (25) ,H2H4 (23,25 
COMMON /02PRP/ 

1 02P1 (13 ) , 02T1 ( 16 ) , 02H1 (13,16) ,02S1(13,16) ,02D1(13,16) , 

2 02P2 ( 13 ) , 02T2 ( 17 ) ,02H2(13,17) ,02S2(13,17) ,02D2(13,17) , 

3 02P3 (5) , 02T3 (61) , 02H3 (5,61) , 02S3(5,61), 02D3(5,61) 
COMMON /H20PRP/ 

1 H20P1 (7 ) , H20T1 ( 13 ) ,H20H1(7,13) ,H20S1(7,13) ,H20D1(7,13) 
COMMON /STD/ 

1 HH2REF , H02REF , HWAREF , SH2REF , S02REF , SWAREF , SH2 A , S02 A , 

2 SWAA 

COMMON /TABLE/ 

1 NH2P (4 ) , NH2T ( 4 ) ,N02P(3) ,N02T(3) ,NH20P(1) ,NH20T(1) 
PARAMETER ( JOULE = 778.16, GC = 32.174 ) 

I FUN C= I FUN C+ 1 


, H2S1 

[15, n; 

, H2D1 

[15,11] 

, H2S2 

20,11 

, H2D2 

20 , 11) 

, H2S3 

29,25 

, H2D3 

29 , 25] 

, H2S4 

23,25 

, H2D4 

[23,25) 


K = 1 

DO 10 1=1 . NBRH 
IF ( IBTYPE ( I ) . EQ . 0 ) THEN 
GO TO 10 
ELSE 

REVM ( I ) =X ( K) **2 
K=K+1 
ENDIF 
CONTINUE 

DO 20 1=1 , NNOD 
IF (INTYPE (I) .EQ.O) THEN 
GO TO 20 
ELSE 

IF (I.EQ.INOZN) THEN 
REVT(I)=X(K) **2 
K=K+1 
ELSE 

REVP (I) =X (K) **2 
REVT ( I ) =X ( K+ 1 ) **2 
K=K+2 
ENDIF 
ENDIF 
CONTINUE 

DO 40 1=1. NPC 
IF (IRTYPE(I) .EQ. 0) THEN 
GO TO 40 
ELSE 

REVR ( I ) =X (K) **2 
K=K+1 
ENDIF 
CONTINUE 

DO 50 1=1 , NRPM 
REVS (I)=X(K) **2 
K=K+1 
CONTINUE 

DO 100 1=1, NNOD 

IF (INTYPE (I) .EQ. 0) THEN 
GO TO 100 
ELSE 

P=REVP ( I ) 

T=REVT ( I ) 


IF i 

( MAT (I 

) .GE. 

3 ) 

GO 

TO 

70 

IF i 

( MAT (I 

) .GE. 

2 ) 

GO 

TO 

60 



QOOO 


c 


60 


70 

72 


C 


74 


80 


C 


90 


C 


c 

100 

c 


110 

c 

c 

120 

c 


130 


CALL PROP ( 1, P, T, 0.0, 0.0, H, S, RHO) 
REVD ( I ) =RHO 
REVH ( I ) =H-HH2REF 
GO TO 100 

CALL PROP ( 2, P, T, 0.0, 0.0, H, S, RHO) 
REVD (I ) =RHO 
REVH ( I ) =H-H02REF 
GO TO 100 

IDHG=1 

IF (IHGN(IDHG) .EQ.I) THEN 
IHG-IDHG 
GO TO 74 
ELSE 

IDHG=IDHG+1 

IF (IDHG.GT.NHGNOD) THEN 
GO TO 74 
ELSE 

GO TO 72 
ENDIF 
ENDIF 


WH2=0 . 0 


DO 80 IH2IN=1 , NH2HG 
IBRH=IH2HG (IHG, IH2I 
WH 2 =WH 2 +RE VM ( IBRH) 
CONTINUE 



W02=0 . 0 

DO 90 I02IN=1 , N02HG ( IHG) 
IBRH=I02HG ( IHG , 102 IN) 

WO 2 =WO 2 +RE VM ( I BRH ) 
CONTINUE 


CEFF=TTB (ICEFF ( IHG) ) 

0F=W02/WH2 

CALL PROP ( 4, P, T, OF, CEFF, H, S, RHO) 
REVD ( I ) =RHO 
REVH ( I ) =H 
ENDIF 


CONTINUE 


F=0 . 0 

DO 120 IMC=1 , NMC 

WIMB ( IMC ) =0 . 6 

DO 110 110=1, NMIO( IMC) 

IOD=MDIR (IMC, IIO) 

IBRH=IMBMC (IMC , IIO) 

WIMB (IMC) =WIMB ( IMC) +IOD*REVM ( IBRH) 
CONTINUE 

F=F+ (WIMB ( IMC) /UMVOL( IMC) ) **2 
CONTINUE 

DO 140 IEC=1,NEC 

EIMB ( IEC) =0 . 6 

DO 130 110=1, NEIO(IEC) 

IOD=EDIR ( IEC , IIO) 

INOD=IEN (IEC, IIO) 

IBRH=IMBEC (IEC, IIO) 

ENTH=REVH(INOD) 

W=REVM ( I BRH ) 

EIMB (IEC) =EIMB (IEC) +IOD*W*ENTH 
CONTINUE 


HEAT=f (RefTl , RefT2 , Ref Ml , RefM2 , BoundaryProps) 
EIMB ( IEC) =EIMB ( IEC) +HEAT 

F=F+ (EIMB ( IEC) /UEVOL( IEC) ) **2 
C 

140 CONTINUE 
C 

DO 160 IPC=1 , NPC 
IBRH=IMBPC (IPC) 

RES=REVR(IPC) 



=REVM(IB 
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c 


c 


c 


c 

170 

C 


C 


180 

C 


190 

C 


INOD5=ITPN(I,5) 
INOD6=ITPN(I, 6) 
ICUR1=ICURV (1,1) 
ICUR2=ICURV (1,2) 
ICUR3=ICURV (1,3) 
ICUR4=ICURV (1,4) 
ICUR5=ICURV (1,5) 
ICUR6=ICURV (1,6) 
ISPEED=ITPS ( I ) 

DIA1=TTB ( IDIA1 ) 
DIA2=TTB(IDIA2 
DIA3=TTB | IDIA3 ) 
W1=REVM ( IBRH1 ) 
W2=REVM ( IBRH2 ) 
W3=REVM ( IBRH3 ) 
P1=REVP ( INOD1 ] 
P2=REVP ( INOD2 ) 
P3=REVP ( INOD3 ) 
P4=REVP ( INOD4 ) 
P5=REVP ( INOD5 j 
P6=REVP ( INOD6 ) 
T1=REVT ( INOD1 ) 
T2=REVT ( INOD2 ) 
T3=REVT ( INOD3 ) 
T4=REVT ( INOD4 ) 
T5=REVT ( INOD5 ) 
T6=REVT ( INOD6 ) 
D1=REVD ( INOD1 ) 
D2=REVD ( INOD2 ) 
D3=REVD ( INOD3 ) 
D4=REVD ( INOD4 ) 

D5— REVD ( INOD5 ) 
D6=REVD(INOD6) 
H1=REVH ( INOD1 ) 
H2=REVH ( INOD2 ) 
H3=REVH ( INOD3 ) 
H4=REVH ( INOD4 ) 
H5=REVH ( INOD5 ) 
H6—REVH ( INOD6 ) 
SPEED==REVS (ISPEED) 


ETP C I1=W1*H1+W2 *H3+W3 *H5-W1*H2-W2 *H4-W3 *H6 

CALL TPIMB (I.ICUR1,ICUR2.DIA1 / W1,SPEED,D1,D2,P1,P2,T1,T2,H1, 
1 H2 ,T1TP(I) , T2TP ri) ) 

CALL TPIMB ( I , ICUR3 , I CUR 4 , DIA2 . W2 , SPEED , D3 , D4 , P3 , P4 , T3 , T4 , H3 , 
1 H4,P1TP(I) ,P2TP(I) ) 

CALL TPIMB (I, ICUR5 , ICUR6 , DIA3 . W3 , SPEED, D5 , D6 , P5 , P6 , T5 , T6 , H5 , 
1 H6 , P3TP ( I ) ,P4*P(I) J 


1 

2 

3 

4 

5 

6 

ENDIF 


F + (ETP (I) /UETP (I) ) 
(TITP(I) /UTITP(I) 
(T2TP ( I ) /UT2TP ( I ) 
( P1TP ( I ) /UP1TP ( I ) 
(P2TP(I)/UP2TP(I) 
(P3TP(I)/UP3TP(I 
(P4TP(I)/UP4TP(I) 


**2 + 
**2 + 
**2 + 
**2 + 
**2 + 
**2 + 
**2 


CONTINUE 


WH2=0 . 0 
W02=0 . 0 


DO 180 IH2 IN=1 , NH2HG ( 
IBRH2=IH2HG ( IHGNOZ . IH 
WH2=WH2+REVM ( IBRH2 ) 
CONTINUE 


IHGNOZ) 

2IN) 


DO 190 102 IN=1,N02HG( IHGNOZ) 
IBRH2=I02HG ( IHGNOZ . I02IN) 
W02=W02+REVM ( IBRH2 ) 


CONTINUE 


ATHROT=TTB ( IATHRT) 
OFCOM =TTB (IMIXR) 
WT0TAL=W02 +WH2 
0FCAL=W02 /WH2 



c 

c 


c 

c 


210 


220 


230 


240 


C 


C 


OFIMB=OFCAL-OFCOM 

F=F+ (OFIMB/UOF) **2 

GAMH2N— 1 . 3 
GAM02N— 1 . 3 
GAMH20=1 . 3 
PTOTAL=REVP ( INOZN ) 

TTOTAL=REVT ( INOZN) 

CEFF=TTB ( I CEFF ( IHGNOZ ) ) 

XF = 1.0 / (1.0 + OFCAL) 

XO = 1.0 - XF 

XH2 = XF - XO * 2.0 * CEFF * 2.016 / 31.9988 

XH20 = XO * 2.0 * CEFF * 18.0153 / 31.9988 

X02 = 1.0 - XH2 - XH20 

GAMMA=XH2 *GAMH2N+XH20*GAMH20+X02 *GAM02N 

GASCON= (XH2/2 . 016+XH20/18 . 0153+X02/31 . 9988) *1545 . 3 

TTHROT=TTOTAL/ ( 1+ ( GAMMA- 1 ) / 2 ) 

PTHROT=PTOTAL/ (1+ (GAMMA-lj/2 ) ** (GAMMA/ ( GAMMA- 1) ) 
VTHROT= S QRT ( GAMMA * GAS CON * TTHROT * G C ) 

RHO=144 . 0*PTHROT/ (GAS CON* TTHROT) 
WNOZ=RHO*ATHROT*VTHROT 
WCAL=REVM ( INOZBR) 

WZ IMB=WCAL-WNOZ 


F=F+ (WZIMB/UMVOL (INOZMC) ) **2 

IF (IFUNC.EQ.l) THEN 
DO 210 1=1 , NMC 
WIMB0 ( I ) =WIMB ( I ) 

CONTINUE 
DO 220 1=1, NEC 
EIMBO ( I ) =EIMB ( I ) 

CONTINUE 

DO 230 1=1, NPC 

PIMBOm=PIMB(I) 

CONTINUE 
DO 240 1=1. NTPC 
IF (ITPTY(I) .EQ.l) THEN 
ETPO(I) =ETP ( I ) 

T1TP0 ( I ) =T1TP ( 1 
T2TP0 (I) =T2TP (I 
P1TP0 ( I ) =P1TP ( I 
P2TP0 (I) =P2TP (I 
ELSE 

ETPO(I) =ETP(I) 

T1TP0 ( I ) =T1TP ( I 
T2TP0 (I) =T2TP (I 
P1TP0 ( I ) =P1TP ( I 
P2TP0 (I) =P2TP (I 
P3TP0 (I) =P3TP (I 
P4TP0 (I) =P4TP (I 
ENDIF 
CONTINUE 
WZ IMB0=WZ IMB 
OFIMBO=OFIMB 


ELSE 

ENDIF 

RETURN 

END 


C 

C 

c***************************************************************** 
SUBROUTINE PROP (MAT, PRSI , TMPI , OF, CEFF, ZENTH, ZENTR, ZDENS) 

C 

C PROP - PROPERTY PROGRAM CALCULATING HYDROGEN, 

C OXYGEN, STEAM AND HOT GAS PROPERTIES 

C 

COMMON /H2PRP/ 

* H2P1 ( 15 ) , H2T1 ( 11) , H2H1 ( 15 , 11) ,H2S1(15,11) ,H2D1(15,11) , 

* H2P2 (20 ,H2T2 (11) ,H2H2 (20,11) ,H2S2 (20,11 ,H2D2 (20,11), 

* H2P3 (29 , H2T3 (25) ,H2H3 (29,25) ,H2S3 (29,25 ,H2D3 (29,25) , 

* H2P4 (23) , H2T4 (25) ,H2H4 (23,25) ,H2S4 (23,25) ,H2D4 (23,25) 

COMMON /02PRP/ 

* 02P1 ( 13 ) , 02T1 ( 16 ) , 02H1 ( 13 , 16) ,02S1(13,16) ,02D1(13,16) , 

* 02P2 ( 13 ) ,02T2 17) ,02H2(13,17j ,02S2(13,17) ,02D2(13,17) , 



02P3 (5) , 02T3(61) ,02H3(5,61) , 02S3(5,61), 02D3(5,61) 
COMMON /H20PRP/ 

H20P1 ( / ) ,H20T1 ( 13 ) ,H20H1(7,13) ,H20S1(7,13) ,H20D1(7,13) 
COMMON /TABLE/ 

NH2P ( 4 ) , NH2T ( 4 ) / N02P(3) ,N02T(3) ,NH20P(1) ,NH20T(1) 
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c 


* 38.340,38.745,39.105,39.419,39.683,39.894,40.049,40.144/ 


C 


C 


c 


c 

c 

51 

52 

53 


C 

c 

c 


c 

c 

10 


c 


11 


c 


12 


c 


13 


C 


14 


DATA (SL02 (J) , J=l, 16) / 

* 0.698,0.708,0.717,0.727,0.736,0.746,0.755,0.764,0.772, 

* 0.781,0.790,0.798,0.806,0.815,0.823,0.831/ 


DATA ( SV02 ( J) , J=1 , 16) / 

* 1.273,1.263,1.254,1.245,1.237,1.229,1.221,1.214,1.207, 

* 1.200,1.193,1.187,1.180,1.174,1.168,1.162/ 


DATA (DL02 (J) , J=l, 16) / 

* 71.630,70.941,70.243,69.536,68.818,68.089,67.347,66.593, 

* 65.823,65.037,64.234,63.412,62.567,61.699,60.804,59.880/ 


DATA (DV02(J) ,J=1,16)/ 

* 0.246,0.305,0.374,0.455,0.547,0.653,0.774,0.911,1.065, 

* 1.239,1.433,1.650,1.893,2.162,2.461,2.794/ 


FORMAT (/3X, 'PROP - REQUESTED PRS > ' F7.2.2X, 

* 'AND TMP > ' ,F7.2,2X ,' FOR H2 IS OUT OF RANGE') 
FORMAT (/3X, 'PROP - REQUESTED PRS > ' F7.2.2X, 

* 'AND TMP > ' ,F7.2,2X, 'FOR 02 IS OUT OF RANGE') 
FORMAT (/3X, 'PROP - REQUESTED PRS > ',F7.2.2X, 

* 'AND TMP > ' , F7 . 2 , 2X, ’ FOR STEAM IS OUT OF RANGE') 

** INTERPOLATE RESULTS FROM SINGLE ARRAY ** 


IPRP=0 
NPX1=2 
NPY1=2 
ZENTH=0 . 0 
ZENTR=0 . 0 
ZDENS=0. 0 


GO TO (10,20,30,40) MAT 


IF (TMPI . GT . 30.0. AND . TMPI . LT . 50.0 
IFfTMPI.GT. 70. 0. AND.TMPI.LT. 110.0 
IFfTMPI.GT. 240. 0. AND.TMPI.LT. 720.0 
IF (TMPI. GT. 1400.0. AND.TMPI.LT. 2 000.0 
GO TO (11,12,13,14) IPRP 


IPRP=1 

IPRP=2 

IPRP—3 

IPRP=4 


IF (PRSI . LT . 20.0. OR . PRS I . GT . 370.0) GO TO 

CALL PRPSAT ( PRSI , TMPI , ZENTH , 

* TSH2 (11) ,NH2P(1) ,NH2T(1) ,11,29.95,50.05, 

* H2P1 , H2T1 , H2H1 , PSH2 , TSH2 , HLH2 , HVH2 ) 

CALL PRPSAT ( PRSI , TMPI , ZENTR , 

* TSH2 (11) ,NH2P(1) , NH2T(1) ,11,29.95,50.05, 

* H2P1 , H2T1 , H2S1 , PSH2 , TSH2 , SLH2 , SVH2 ) 

CALL PRPSAT (PRSI , TMPI , ZDENS , 

* TSH2 (11) , NH2P(1) ,NH2T(1) ,11,29.95,50.05, 

* H2P1 , H2T1 , H2D1 , PSH2 , TSH2 , DLH2 , DVH2 ) 
RETURN 


50 


IFfPRSI.LT. 3400.0. OR. PRSI. GT. 7200.0) GO TO 50 
CALL ITERP2 ( PRSI , TMPI , H2P2 , H2T2 , H2H2 , 

* NH2 P ( 2 ) , NH2T ( 2 ) , NPX1 , NPY1 , NH2 P f 2 ) . ZENTH , N1 ) 
CALL ITERP2 (PRSI , TMPI , H2P2 , H2T2 , H2S2 , 

* NH2Pf2] .NH2T(2) , NPX1 , NPY1 , NH2P ( 2) , ZENTR, Nl) 
CALL ITERP2 (PRSI , TMPI , H2P2 , H2T2 , H2D2 , 

* NH2P ( 2 ) ,NH2T(2) , NPX1 , NPY1 , NH2P ( 2 ) , ZDENS, Nl) 
RETURN 


IFfPRSI.LT. 1400.0. OR. PRSI. GT. 7000.0) GO TO 50 
CALL ITERP2 ( PRS I , TMPI , H2 P3 , H2 T3 , H2H3 , 

* NH2 P ( 3 ) , NH2T ( 3 ) , NPX1 , NPY 1 , NH2 P ( 3 ) , ZENTH , Nl ) 
CALL ITERP2 (PRSI. TMPI ,H2P3 ,H2T3 .H2S3 , 

* NH2P( 3 ) , NH2T ( 3 ) , NPX1 , NPY 1 , NH2 P (3 ) , ZENTR, Nl) 
CALL ITERP2 ( PRS I , TMPI ,H2P3,H2T3.H2D3. 

* NH2P ( 3 ) , NH2T ( 3 ) , NPX1 , NPY 1 , NH2 P ( 3 ) , ZDENS, Nl) 
RETURN 

IFfPRSI.LT. 1400.0. OR. PRSI. GT. 5800.0) GO TO 50 
CALL ITERP2 (PRSI , TMPI , H2P4 , H2T4.H2H4, 

* NH2P(4 ) ,NH2T(4) , NPX1 , NPY 1 , NH2 P (4 ) , ZENTH, Nl) 
CALL ITERP2 (PRS I .TMPI .H2P4 ,H2T4 ,H2S4 , 

* NH2P ( 4 ) , NH2T ( 4 ) , NPX1 , NPY1 , NH2P ( 4 ) , ZENTR, Nl) 



ooo ooo o on no noon on no o on 


CALL ITERP2 (PRSI,TMPI,H2P4 ,H2T4 ,H2D4 , 

* NH2P ( 4 ) , NH2T (4 ) , NPX1 , NPY1 , NH2P ( 4 ) ,ZDENS,N1) 
RETURN 


20 IF (TMPI . GT . 160. 0. AND.TMPI.LT. 240.0) IPRP=1 
IF (IPRP. EQ. 1. AND. PRSI. LT. 650.0) IPRP=1 
IF(IPRP.EQ. 1. AND. PRSI.GT. 650.0) IPRP=2 

IF (TMPI. GT. 600.0. AND.TMPI.LT. 1500.0) IPRP=3 
GO TO (21, 22 , 23) IPRP 

21 IF (PRSI . LT . 30.0. OR . PRS I . GT . 630.0) GO TO 50 

IF (TMPI. LT. 160.0. OR. TMPI. GT. 219.9) GO TO 50 
CALL PRPSAT (PRSI. TMPI, ZENTH, 

* TS02 ( 16) ,N02P(1) ,N02T(1) ,16,159.95.220.05, 

* 02P1 , 02T1 , 02H1 , PS02 , TS02 , HL02 , HV02 ) 

CALL PRPSAT (PRSI. TMPI. ZENTR. 

* TS02 (16) , N02P(iJ ,N02T(1) ,16,159.95,220.05, 

* 02P1 , 02T1 , 02S1 , PS02 , TS02 , SL02 , SV02 ) 

CALL PRPSAT (PRS I , TMPI . Z DENS . 

* TS02 (16) , N02P ( i) , N02T(1) ,16,159.95,220.05, 

* 02P1 , 02Tl , 02D1 , PS02 , TS02 , DL02 , DV02 ) 

RETURN 

22 IF(PRSI.LT. 2000.0. OR. PRSI. GT. 8000.0) GO TO 50 
CALL ITERP2 ( PRSI , TMPI , 02P2 , 02T2 , 02H2 , 

* N02P(2) ,N02T(2) ,NPXl , NPYl , N02P (2 ) , ZENTH, Nl) 
CALL ITERP2 ( PRSI . TMPI , 02P2 , 02T2 , 02S2 , 

* N02P (2 ) , N02T (2 ) . NPX1 , NPYl , N02P ( 2 ) , ZENTR, Nl) 
CALL ITERP2 (PRSI , TMPI , 02P2 , 02T2 , 02D2 , 

* N02P (2 ) , N02T ( 2 ) , NPX1 , NPYl , N02P ( 2 ) ,ZDENS,N1) 
RETURN 

23 IF(PRSI.LT. 2000.0. OR. PRSI. GT. 4000.0) GO TO 50 
CALL ITERP2 (PRSI , TMPI , 02P3 , 02T3 , 02H3 , 

* N02PX3) ,N02T(3) ,NPX1,NPY1,N02P(3) , ZENTH, Nl) 
CALL ITERP2 (PRSI, TMPI, 02P3 , 02T3 , 02S3 , 

* N02PX3).N02T(3) , NPX1 , NPY 1 , N02 P ( 3 ) , ZENTR , Nl ) 
CALL ITERP2 (PRSI , TMPI , 02 P3 , 02T3 .02D3 , 

* N02P ( 3 ) , N02T ( 3 ) , NPXl , NPYl , N02P ( 3 ) ,ZDENS,N1) 
RETURN 


30 IF(TMPI.LT. 1400.0. OR. TMPI. GT. 2000.0) GO TO 50 
IF (PRSI . LT . 100.0. OR. PRSI. GT. 700.0) GO TO 50 
CALL ITERP2 ( PRSI , TMPI . H20P1, H20T1 , H20H1 , 

* NH20P ( 1) , NH20T ( 1) , NPXl, NPYl, NH20P ( 1) , ZENTH, Nl) 
CALL ITERP2 (PRSI , TMPI.H20P1 , H20T1 , H20S1 , 

* NH20P ( 1) , NH20T ( 1) , NPXl , NPYl , NH20P ( 1] , ZENTR, Nl) 
CALL ITERP2 (PRSI , TMPI ,H20P1, H20T1 ,H20 d1, 

* NH20P ( 1) , NH20T ( 1) , NPXl , NPYl , NH20P ( 1) ,ZDENS,N1) 
RETURN 


4 0 CALL PRPMIX ( PRSI , TMPI , OF , CEFF , HMIX , SMIX) 
************ MODIFIED 2/19/93 ********** 

4 0 CALL PRPMIX ( PRSI , TMPI , OF , CEFF , HMIX , SMIX , DMIX) 
ZENTH=HMIX 
ZENTR=SMIX 
ZDENS=0 . 0 

************ MODIFIED 2/19/93 ************ 
ZDENS=DMIX 
RETURN 


50 


IF (MAT. EQ. 1) 
IF (MAT. EQ. 2) 
IF (MAT. EQ. 3) 
RETURN 


WRITE (21,51) 
WRITE (21,52) 
WRITE (21, 53) 


PRS I, TMPI 
PRS I, TMPI 
PRSI , TMPI 


END 

**************************************************************** 
SUBROUTINE PRPMIX (P, TMPI , OF, CEFF, HMIX, SMIX) 

************ MODIFIED 2/19/93 ************ 

SUBROUTINE PRPMIX (P, TMPI , OF , CEFF, HMIX, SMIX, DMIX) 

PRPMIX - CALCULATES HOT GAS MIXTURE PROPERTIES. 



c 


c 

c 


c 


c 


c 


c 

c 

c 


c 


c 


COMMON /H2PRP/ 

* H2P1 (15) , H2T1 (11) ,H2H1 (15 , 11) ,H2S1(15,11) ,H2D1(15,11) , 

* H2P2 (20) ,H2T2 ill) ,H2H2(20,11) ,H2S2(20,11) ,H2D2(20,11) , 

* H2P3 (29 ) ,H2T3(25 ,H2H3 (29,25) ,H2S3(29,25 ,H2D3 (29,25) , 

* H2P4 (23 ) , H2T4 (25) ,H2H4(23,25) ,H2S4 (23,25) ,H2D4 (23,25) 
COMMON /02PRP/ 

* 02P1 ( 13 ) , 02T1 (16) , 02H1 (13,16) ,02S1(13,16) ,02D1(13,16) , 

* 02P2 (13 ) , 02T2 ( 17 ) ,02H2(13,17) ,02S2(13,17) ,02D2(13,17) , 

* 02P3 (5) , 02T3 (61) , 02H3 (5,61) , 02S3(5,61), 02D3(5,61) 
COMMON /H20PRP/ 

* H20P1 (7 ) , H20T1 (13 ) ,H20H1(7,13) ,H20S1(7,13) ,H20D1(7,13) 
COMMON /TABLE/ 

* NH2P (4 ) ,NH2T (4 ) ,N02P(3) ,N02T(3) ,NH20P(1) ,NH20T(1) 
COMMON /STD/ 

* HH2REF , H02REF , HWAREF , SH2REF , S02REF , SWAREF , SH2 A , S02A , 

* SWAA 


XMWH2 = 2.0160 

XMW02 = 31.9988 

XMWH20 = 18.0153 

HCOMB = -6825.6550 

NPX1 = 2 
NPY1 = 2 
ITST1= 0 
ITST2= 0 
ITST3= 0 
ITST4= 0 
ITST5= 0 
ITST6= 0 

XF = 1.0 / (1.0 + OF) 

XO = 1.0 - XF 

XH2 = XF - XO *2.0* CEFF * XMWH2 / XMW02 
XH20 = XO * 2.0 * CEFF * XMWH20 / XMW02 

X02 = 1.0 - XH2 - XH20 

EH2 = XH2 / XMWH2 
EH20 = XH20 / XMWH20 
E02 = X02 / XMW02 

ET = EH2 + EH20 + E02 

YH2 = EH2 / ET 
YH20 = EH20 / ET 
Y02 = 1.0 - YH2 - YH20 

PH 2 = P * YH2 

PH20 = P * YH20 

P02 = P * Y02 

IF(TMPI.LT. 1000.0. OR. TMPI.GT. 2000.0) ITST1=1 
IFjPH2 .LT. 1400.0. OR. PH2 .GT. 5800.0) ITST2=1 
CALL ITERP2 (PH2 , TMPI , H2 P4 , H2T4 , H2H4 , 

* NH2P(4) ,NH2T(41,NPX1,NPY1,NH2P(4) ,HH2,N1) 
CALL ITERP2 ( PH 2 , TMPI , H2 P4 , H2T4 , H2S4 , 

* NH2P ( 4 ) , NH2T(4) , NPX1 , NPY1 , NH2P (4 ) ,SH2,N1) 

IF (TMPI. LT. 1400. 0. OR. TMPI.GT. 2000.0) ITST3=1 
IF (PH20. LT . 100.0. OR. PH20.GT. 700.0) ITST4=1 
CALL ITERP2 fPH20,TMPI .H20P1 .H20T1,H20H1, 

* NH20P ( lj , NH20T ( l) , NPX1 , NPY1 , NH20P (11 , HH20 , N1 ) 
CALL ITERP2 ( PH20 , TMPI , H20P1 , H20T1 , H20S1 , 

* NH20P ( 1) , NH20T ( 1) , NPX1 , NPY1 , NH20P ( 1) ,SH20,N1) 

IF (Y02.LT. 0.001) THEN 
DH02 =0.0 
DS02 =0.0 
ELSE 

IF (TMPI . GT . 600.0. AND. TMPI . LT . 1500 . 0 ) ITST5=1 
IF]P02 . LT .2000.0. OR . P02 .GT. 4000.0) ITST6=1 
CALL ITERP2 (P02 , TMPI , 02P3 , 02T3 , 02H3 , 

* N02P(3) ,N02T(3) , NPX1 , NPY1 ,N02P ( 3 ) ,H02,N1) 

CALL ITERP2 (P02 , TMPI , 02 P3 , 02 T 3 , 02 S3 , 

* N02P ( 3) ,N02T(3) ,NPX1,NPY1,N02P(3) ,S02,N1) 

DH02 = HO 2 - H02REF 

DS02 = S02 - S02REF + S02A 



ENDIF 

C 

10 DHH2 = HH2 - HH2REF 

DHH20M = (HH20 - HWAREF) + HCOMB 
DSH2 = SH2 - SH2REF + SH2A 
DSH20 = SH20 - SWAREF + SWAA 
C 

HMIX = XH2*DHH2 + XH20*DHH20M + X02*DH02 

SMIX = XH2*DSH2 + XH20*DSH20 + X02*DS02 

C 

C ************ MODIFIED 2/19/93 ************ 

XMWMIX = XMWH20*YH20 + XMWH2*YH2 + XMW02*Y02 
DMIX = 144 . *XMWMIX*P/ ( 1545 . 3*TMPI) 

C 

IF (ITST1.EQ. l.OR. ITST2 . EQ. 1) WRITE(21,51) PH2.TMPI 
IF (ITST3.EQ.1.0R.ITST4.EQ.1 WRITE(21,52 PH20,TMPI 
IF (ITST5.EQ. l.OR. ITST6.EQ. 1) WRITE(21,53) P02,TMPI 
C 

51 FORMAT (/3X, 'PRPMIX - REQUESTED PH2 PRS > * ,F7.2,2X, 

* 'AND TMP > ' F7.2,2X , ' FOR " H2" IS OUT OF RANGE') 

52 FORMAT (/3X, 'PRPMIX - REQUESTED PH20 PRS > ' ,F7.2,2X, 

* 'AND TMP > ' F7.2,2X ,'FOR "H20" IS OUT OF RANGE') 

53 FORMAT (/ 3 X , • PRPMIX - REQUESTED P02 PRS > ' ,F7.2,2X, 

* 'AND TMP > ' ,F7.2,2X, 'FOR " 02" IS OUT OF RANGE') 

C 

RETURN 

END 

C** ************** ************************************** *********** 


c 

c 

c 

c 


SUBROUTINE PRPSAT 
* PRS 1 , TMP1 , PROP , P: 


X , Y , FPROP . TCRT , NX1 , NY1 , NX2 , YL, YH 
:S2 , TMP2 , PROPL, PROPV) 


t 


PRPSAT - CALCULATES NBS PROPERTIES NEAR SATURATION CURVE 


DIMENSION PRS 1(1), TMP1 ( 1 ) 

NR1=NX1 

NPX1=2 

NPY1=2 

NPX2=2 


ZPLGAS-O. 0 
ZPHGAS=0 . 0 
ZPLLIQ=0.0 
ZPHLIQ=0 . 0 
ZPROP1=0 . 0 
ZPROP=0 . 0 
FPROP=0 . 0 
ZTSAT=0 . 0 
ARGA=0 . 0 


ARGB=0 . 0 
ZTSATT=0 . 0 


C 

c 

CALL ITERP2 (X, Y , PRS1 , TMP1 , PROP , NX1 , NY1 , NPX1 , NPY1 , NR1 , ZPR0P1 , N1 ) 

FPR0P=ZPR0P1 

IF f Y . GT . TCRT) GO TO 70 

CALL ITERP1 (X, PRS 2 , TMP 2 , NX 2 , NPX2 , ZTSAT , N2 ) 

IF(Y.LT.ZTSAT) GO TO 61 
C 

C * * GAS CALCULATIONS * * 

C 

CALL ITERP1 ( X , PRS 2 , PROPV . NX2 , NPX2 , Z PGAS , N2 ) 

CALL ITERP2 (X . YH, PRS1 , TMP1 , PROP, NX1 , NY1, NPX1 , NPY1 , NR1 , ZTST, Nl) 

DTST=ZTST-ZPGAS 

IFfDTST.GT. 0.0001) GO TO 50 

ZPLGAS=ZPGAS 

IF (ZPROP1 . LT. ZPGAS) GO TO 70 
GO TO 51 

50 Z PHGAS= Z PGAS 

IF (ZPROP1.GT. ZPGAS) GO TO 70 
C 

51 LPR=1 

53 PRSD=PRS1 (LPR) -0.0001 
IF(PRSD.GT.X) GO TO 52 
LPK— LPR+1 
GO TO 53 


C 

52 ARGA=PRS1 (LPR) 
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c 


CALL ITERP1 (ARGA, PRS2 , TMP2 , NX2 , NPX2 , ZTSATT , N2 ) 


C 


C 


LTP=1 

54 TMPD=TMP1 (LTP) -0. 0001 

IF fTMPD.GT. ZTSATT) GO TO 55 

LTP=LTP+1 

GO TO 54 


55 ARGB=TMP1 (LTP) 

YY=ARGB 

IF(DTST.GT. 0.0001) CALL ITERP2 ( X , Y Y , PRS 1 , TMP1 , PROP , NX1 , NY 1 , 

* NPX1 . NPY1 . NR1 . ZPLGAS , Nl) 

IF (DTST . LT .6.0001) CALL ITERP2 ( X , YY , PRS 1 , TMP1 , PROP , NX1 , NY 1 , 

* NPX1,NPY1.NR1,ZPHGAS,N1) 

ZPROP=ZPHGAS- (ZPHGAS- ZPLGAS) * ( (ARGB-Y) / (ARGB-ZTSAT) ) 

FPROP=Z PROP 


GO TO 70 


* * LIQ CALCULATIONS * * 

61 CALL ITERP1 (X , PRS2 , PROPL, NX2 , NPX2 , ZPLIQ , N2) 

CALL ITERP2 ( X , YL , PRS 1 , TMP1 , PROP , NX1 , NY1 , NPX1 , NPY 1 , NR1 , ZTST , N1 ) 
DTST=ZTST- ZPLIQ 
IF(DTST.GT. 0.0001) GO TO 59 
ZPLLIQ=ZPLIQ 

IF (ZPROP1.LT. ZPLIQ) GO TO 70 
GO TO 60 

59 Z PHLI Q= Z PLI Q 
IF(ZPR0P1.GT. ZPLIQ) GO TO 70 

60 LPR=1 

63 PRSD=PRS1(LPR) -0.0001 
IF(PRSD.GT.X) GO TO 62 
LPR=LPR+ 1 

GO TO 63 

62 ARGA=PRS 1 ( LPR- 1 ) 

CALL ITERP1 ( ARGA , PRS 2 , TMP2 , NX2 , NPX2 , ZTSATT , N2 ) 

LTP=1 

64 TMPD=TMP1(LTP) -0.0001 

IF fTMPD.GT. ZTSATT) GO TO 65 

LTP=LTP+1 

GO TO 64 

65 ARGB=TMP 1 ( LTP- 1 ) 

YY=ARGB 

IF (DTST.GT. 0.0001) CALL ITERP2 ( X , YY , PRS 1 , TMP1 , PROP , NX1 , NY 1 , 

* NPX1 , NPY1 , NR1 , ZPLLIQ , Nl) 

IFfDTST.LT. 6.0001) CALL ITERP2 (X, YY , PRS1 , TMP1 , PROP, NX1 , NY1 , 

* NPX1 , NPY1 , NR1 , ZPHLIQ , Nl) 

ZPROP=ZPHLIQ- ( ZPHLIQ-ZPLLlQ) * ( (ZTSAT-Y) / (ZTSAT-ARGB) ) 

F PRO P= Z PRO P 


70 CONTINUE 


RETURN 

END 

**************************************************************** 
SUBROUTINE ITERP1 ( X , XT , YT , NX , NPX , Y , NERR) 

ITERP1 - SINGLE INTERPOLATION ROUTINE. 

DIMENSION XT ( 1 ) , YT ( 1 ) 

NERR-0 

INTER=1 

NP=NPX 

IF (NX .LT. NP) NP=NX 
IH=NP/2 
1=1 

IF (XT (I) -X) 30,20,10 
10 IH=0 

12 NERR=1 
GO TO 70 

13 NERR=2 
GO TO 70 

20 INTER=2 



22 Y=YT(I) 

GO TO 999 
30 I=NX 

IF(XT(I)-X) 13,20,40 
40 Nl=l 
N2=NX 

45 MP= (N1+N2 ) /2 

50 IF (XT (MP) -X) 52,54, 56 

52 N1=MP 

GO TO 60 
54 I=MP 

GO TO 20 
56 N2=MP 

60 IF ( (N2-N1) .NE. 1) GO TO 45 
C 

IF (N2.GT. (IH+1) ) GO TO 65 
I=IH+1 
GO TO 70 
65 I=N2 
C 

IF (N2 .GT. I) I=N2 
70 K=I-IH 
N=K+NP-1 
Y=0. 

IF(N-NX) 90,90,80 
80 N=NX 

K=NX-NP+1 
90 DO 120 J=K, N 
P=1.0 

DO 110 I=K, N 
IF (I“J) 100,110,100 
100 P=P* rx-XT(I) )/ (XT(J) -XT(I) ) 

110 CONTINUE 

Y=Y+YT (J) *P 
120 CONTINUE 
GO TO 999 

ENTRY ENTERP (X,XT,YT,Y) 

Y=0 . 

GO TO (90,22) , INTER 
999 CONTINUE 
RETURN 
END 

C* **************************************************************** 
SUBROUTINE ITERP2 (X, Y , XT , YT, ZT, NX, NY , NPX, NPY ,NR, Z , NERR) 

C 

C ITERP2 - DOUBLE INTERPOLATION ROUTINE. 

C 

DIMENSION XT ( 1) , YT ( 1) ,ZT(NR,1) ,ZC(15) 

NERRB=0 

NPYY=NPY 

IF (NY .LT. NPY) NPYY=NY 
IH=NPYY/2 
1=1 

IF ( YT (I) -Y) 30,20,10 
10 IH=0 

12 NERRB=201 
GO TO 70 

13 NERRB=204 
GO TO 70 

20 CALL ITERP1 ( X , XT , ZT ( 1 , I ) , NX , NPX , Z , NERRA) 

GO TO 999 
30 I=NY 

IF (YT (I) -Y) 13,20,40 
40 Nl=l 
N2=NY 

45 MP= (N1+N2 ) /2 
50 IF ( YT (MP) -Y) 52,54,56 
52 N1=MP 
GO TO 60 
54 I=MP 

GO TO 20 
56 N2=MP 

60 IF ( (N2-N1) .NE. 1) GO TO 45 
I=N2 

IFfI .LT. (IH+1)) I=IH+1 
70 K=I-IH 

N=K+NPYY-”1 
IF (N-NY) 90,90,80 
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80 N=NY 

K=NY-NPYY+1 
90 J=0 

DO 100 I=K, N 
J=J+1 

IF f J .NE. 1) GO TO 95 

CALL ITERP1(X,XT,ZT(1,I) ,NX,NPX, ZC(J) , NERRA) 
GO TO 100 

95 CALL ENTERP (X, XT, ZT(1,I) , ZC ( J) ) 

100 CONTINUE 

CALL ITERP1 ( Y , YT ( K) , ZC,NPYY,NPYY, Z ,NERRC) 

999 NERR=NERRA+NERRB 
RETURN 
END 


SUBROUTINE OPT (N, NMAX, EPSC , EPST.RHO, ALPHA, BETA, 

1 ISTAGE , F2 , GNORM, DXNORM, X , X2) 

DIMENSION X ( 2 5 ) , X2 (25) ,G(25) , G2 (25), DX (25) ,H(25,25) 

ISTAGE=STAGE COUNTER 

ISTAGE=1 

IFUNC=0 

CALL OBJF (IFUNC, X, F) 

CALL GRAD (IFUNC , N, X, G, F) 

SET INITIAL METRIC H TO THE IDENTITY MATRIX 
10 DO 12 1=1, N 
DO 12 J=1,N 
IF (J.NE.I) THEN 
H (I , J) =0 . 0 

ELSE 

H (I , J) =1 . 0 

ENDIF 

12 CONTINUE 

KC=STEP DIRECTION PARAMETER 

KC=1 NEGATIVE GRADIENT STEP DIRECTION 
KC=0 QUASI -NEWTON STEP DIRECTION 


KC=1 


15 


18 

20 


CALCULATE SEARCH DIRECTION VECTOR DX=-H*G 
DO 20 1=1, N 
DX (I) =0 . 0 
DO 18 J=1 , N 

DX (I ) =DX ( I ) -H ( I , J ) *G ( J) 

CONTINUE 


CALL UNIVARIATE SEARCH ROUTIVE 
70 CALL ARMIJO (IFUNC, I FLAG, N,RHO, ALPHA, BETA, DMAX, 
1 DX , X , X2 , F , F2 . G) 

CALL GRAD (IFUNC, N, X2 ,G2 , F2) 

DO 50 1=1, N 

50 DXDG=DXDG+(X2 (I) -X(I) ) * (G2 (I) -G(I) ) 

IF (DXDG .LT. 0.0 .AND. KC .EQ. 0) IFLAG=1 
IF (IFLAG .NE. 1) GO TO 90 
76 IF (KC .EQ. 1 .AND. IFLAG .EQ. 1) GO TO 200 
GO TO 10 

GRADIENT NORM* *2 < EPSC IMPLIES CONVERGENCE 
90 GNORM=0 . 0 

DO 92 1=1, N 

9 2 GNORM=GNORM+G2 ( I ) * * 2 
IF (GNORM-EPSC) 200,110,110 

DELTA X NORM* *2 < EPST IMPLIES TERMINATION 

10 DXNORM=0 . 0 
DO 112 1=1, N 

12 DXNORM=DXNORM+ ( X2 ( I ) -X ( I ) ) * * 2 
IF (DXNORM-EPST) 200,120,120 

INCREMENT STAGE COUNTER 
120 ISTAGE=ISTAGE+1 
C 



c 

c 

c 


c 

c 


125 

c 

c 

200 

c 

c 

c 


c 

110 

c 


111 

c 


c 

112 

115 

C 


120 


180 

200 

C 

c 


c 

c 

c 


c 

c 

10 

310 


320 

C 

400 

c 


CHECK THAT MAXIMUM NUMBER OF STAGES NOT EXCEEDED 
IF (ISTAGE .GT. NMAX) GO TO 200 

UPDATE METRIC H 

CALL BFGS (IFLAG , N, X, X2 , G, G2 , H) 

IF (IFLAG. EQ.l) GO TO 76 
KC=0 


REINITIALIZE 

F=F2 

DO 125 1=1 ,N 


GO TO 15 


RETURN 


END 


SUBROUTINE ARMIJO (IFUNC. IFLAG, N.RHO, ALPHA, BETA, DMAX, 
1 DX. X, X2 , F. F2 , G) 

DIMENSION X(25) ,X2 (25) ,DX(25) ,G(25) 

DMAX=0 . 05 
IFLAG=0 
ICOUNT=l 
GDX=0 . 0 

DO 110 1=1, N 
GDX=GDX+G (I) *DX (I) 

RAT I 0=0 . 0 
DO 111 1=1, N 
RATI02=ABSJDX (I) /X (I) ) 

RATIO=MAX (RATIO , RATI02 ) 

RMU=RHO 

SCALE=RMU*RATIO 

IF (SCALE .GT. DMAX) RMU=DMAX/RATIO 

DO 115 1=1, N 

X2 ( I ) =X ( I ) +RMU*DX ( I ) 

CALL OBJF ( IFUNC, X2,F2) 

T B AR=F2 - F -RMU * ALPHA * G DX 

IF (TBAR .GT. 0.0) GO TO 120 

IF (F-F2) 180,180,200 

RMU=RMU*BETA 

I COUNT=I COUNT+ 1 

IF (ICOUNT .LE. 12) GO TO 112 

IFLAG=1 

RETURN 

END 


SUBROUTINE BFGS ( IFLAG, N, X, X2 ,G,G2,H) 

DIMENSION X (25 ) ,X2(25) ,G(25) ,G2(25) ,DX(25) ,DG(25) ,HDG(25) , 
1 V (25) ,H (25, 25) 

IFLAG=0 

CALL DFP METRIC UPDATE 

CALL DFP ( N , DXDG , DGHDG , X , X2 , G , G2 , DX , DG , HDG , H) 

IF (DGHDG .GT. 0.0) GO TO 10 

IFLAG=1 

GO TO 400 

DETERMINE BFGS METRIC UPDATE 
DO 310 1=1, N 

V ( I ) = DGHDG **0.5* (DX (I) / DXDG-HDG (I) / DGHDG ) 

DO 320 1=1, N 
DO 320 J=l. N 
H ( I , J ) =V ( I ) * V ( J ) +H ( I , J ) 

RETURN 

END 
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100 


400 

410 


420 


SUBROUTINE DFP (N , DXDG, DGHDG , X , X2 , G, G2 , DX, DG.HDG, H) 
DIMENSION X(25) ,X2 (25},gJ25) ,G2 (25) ,DX(25) ,DG(25) ,600(25) , 
1 DGH (2d) ,11(25,25) 

DO ~ “ " 


(I) =X2 

- 1 . » 1 1 

;t j -x(!) 

Cl) =621 

I)-G(I) 


(I)=0.0 
(I)=0.0 
400 J=1 , N 


DGHDG=0 . 0 
DO 410 1=1, N 
HDG (I) =0 . 0 
DGH “ 

DO 

HDG ( I ) =HDG (I) +H (I , J) *DG ( J) 
DGH ( I j =DGHj I +DG( J ) *H ( J , I ) 
DXDG=DXDG+DX ( I) *DG ( I ) 
DGHDG=DGHDG+DGH (I) *DG(I) 

DO 420 1=1, N 
DO 420 J=1,N 

H (I , J* 

RET 
END 


=H ( I , J ) + DX ( I ) * DX ( J ) / DXDG-HDG ( I ) * DGH (J) / DGHDG 


10 


20 


SUBROUTINE GRAD (IFUNC,N, X, G, F) 
DIMENSION X(25) ,X2 (25) ,G(25) 
DX=0 .001 


DO 10 1=1, N 
X2(I)=X(I) 

DO 20 1=1, N 

X2 (I) =X(I) +DX 

CALL OBJF ( IFUNC , X2 , F2 ) 

G(I)=(F2-F)/DX 

X2(I)=X(I) 

RETURN 

END 


SUBROUTINE TPIMB ( I , ICURA, ICURB, DIA, W, SPEED, D1 , D2 , PI , P2 , T1 , T2 , 

1 Hi , H2 , ZIMB1 , ZIMB2 ) 

***** TURBINE CURVES REQUIRE ICURA & ICURB < 20 ******************* 

IF XICURA.LT.20) THEN 
F LOWC=W * S QRT (Tl) /PI 
TMACH=SPEED/SQRT (Tl) 

CALL CURVES ( ICURA , FLOWC , TMACH , CHR3 , CHR4 , PRATIO ) 

CALL CURVES ] ICURB , FLOWC , TMACH , CHR3 , CHR4 , TRATIO ) 
DPCHR=P1-P1/PRATI0 
DPCAL=P1-P2 
Z IMB1=DPCAL-DPCHR 
DTCHR=Tl-ri. O-TRATIO) *T1 
DTCAL=T1-T2 
ZIMB2=DTCAL-DTCHR 
ELSE 

***** PUMP CURVES REQUIRE ICURA & ICURB >= 20 ********************* 
FLOWC=W/ (D1*SPEED*DIA**3 ) 

CALL CURVES ( ICURA , FLOWC , CHR2 , CHR3 , CHR4 , HEADC ) 

CALL CURVES j ICURB , FLOWC , CHR2 , CHR3 , CHR4 , EFFCHR) 

DPCHR=HE ADC *D1*SPEED**2*DIA**2 
DPCAL=P2 -PI 
Z IMB1=DPCAL-DPCHR 
PWR=W * (H2-H1) 

PWRID= (144. 0/778. 16) * (W* (P2-P1) /Dl) 

EFFCAL=PWRID/PWR 
ZIMB2= (EFFCAL-EFFCHR) *PWR 
ENDIF 
C 

10 CONTINUE 
C 

RETURN 

END 



